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The  objective  of  this  research  program  is  to  develop  an  optimum 
multi-discriminant/detection  procedure  for  earthquake  and  under- 
ground explosions  with  emphasis  on  events  occurring  within  the 
Asian  continent.  The  approach  to  the  seismic  discrimination 
problem  incorporates  a number  of  diverse  topics  including:  ex- 

plosion and  earthquake  source  modeling?  stress  wave  propagation 
through  realistic  earth  structures  and  prediction  of  teleseis- 
mically  recorded  ground  motion;  the  development  of  signal \ 
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FOREWORD 

This  semi-annual  technical  report  entitled,  "Applica- 
tion of  Advanced  Methods  for  Identification  and  Detection 
of  Nuclear  Explosions  from  the  Asian  Continent,  is  sub 
mitted  by  Systems,  Science  and  Software  (S3)  to  the  Advanced 
Research  Projects  Agency  and  to  the  Air  Force  Office  of 
Scientific  Research  (AFOSR) . This  report  presents  the  results 
of  a continuing  effort  to  obtain  an  optimum  multi-discriminant/ 
detection  procedure  for  earthquakes  and  explosions  occurring 
within  the  Asian  Continent.  The  work  is  being  performed  under 
Contract  Number  F44620-74-C-0063 . Mr.  William  J.  Best  is  the 
AFOSR  technical  contracting  officer. 

Dr.  J.  Theodore  Cherry  is  the  S3  project  manager.  Drs. 
Thomas  C.  Bache  and  Joseph  F.  Masao  are  responsible  for  the 
development  and  application  of  the  seismic  ground  motion  pre- 
diction work.  Dr.  John  M.  Savino  and  MeSBrB.  Kenneth  G._ 
Hamilton  and  David  G.  Lambert  are  responsible  for  the 
analysis  of  the  seismic  data  against  which  all  theoretical 
developments  must  eventually  be  tested.  Acting  as  consul- 
tants on  the  project  are  Professors  Charles  B.  Archambeau 
of  the  University  of  Colorado,  David  G.  Harkrider  of  the 
California  Institute  of  Technology  and  Donald  V.  Helmberger 
of  the  California  Institute  of  Technology. 

The  authors  wish  to  extend  their  sincere  appreciation 
to  Ms.  Bernadine  Ludwig  and  Ms.  Darlene  Roddy  for  the  many 
hours  spent  on  the  preparation  of  this  report. 
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I .  INTRODUCTION 

The  conduct  of  the  research  under  this  contract  is 
aimed  at  the  development  of  an  optimum  multi-discriminant/ 
detection  procedure  for  earthquakes  and  underground  explo- 
sions with  emphasis  on  events  occurring  within  the  Asian 
continent.  In  order  to  realize  this  objective  our  research 
program  has  involved  a combined  theoretical/empirical  ap- 
proach. That  is,  deterministic  predictions  of  teleseismic 
ground  motion  generated  by  underground  explosion  and  earth- 
quake sources  are  ultimately  compared  to  actual  observations. 
In  this  way  we  are  provided  with  a confirmed  theoretical 
framework  for  testing  existing  discriminants,  as  well  as 
designing  new  discriminants. 

Our  approach  is  quite  comprehensive  and,  in  outline 
form,  involves  the  following: 

1.  Explosion  and  earthquake  source  modeling. 

2.  Stress  wave  propagation  through  complicated 
realistic  earth  structure. 

3.  Development  of  state  of  the  art  signal  enhance- 
ment and  identification  techniques. 

4.  Multi-discriminant  design  and  evaluation. 

Section  II  of  this  report  is  devoted  to  a detailed 
description  of  an  extremely  versatile  computer  code,  MARS 
(Multiple  Arrival  Recognition  System) , that  incorporates 
several  different  data  processing  techniques  for  signal  de- 
tection, enhancement  and  identification.  The  MARS  code  now 
accepts  from  one  up  to  three  components  of  seismic  data  and 
applies  a series  of  narrow  band  filters  to  determine  spectral 
amplitudes.  Corrections  for  instrument  response  may  be 
made  and  true  dispersion  data  computed.  Polarization  and 
dispersion  filters  may  then  be  applied  to  separate 
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the  various  modes  of  wave  propagation  present  in  the  record. 
Then,  for  example,  the  code  can  filter  the  record  to  select 
only  those  P-waves  arriving  along  selected  ray  paths  (azi- 
muth and  emergence  angle).  Matched  filtering  and  cross- 
correlation are  among  the  other  data  processing  capabilities 
available  as  options  in  the  code. 

The  variable  frequency  magnitude  (VFM)  technique, 
embodied  in  the  MARS  code,  has  been  tested  on  a wide  range 
of  seismic  event  data  to  determine  its  effectiveness  as  a 
discriminant  between  earthquakes  and  underground  explosions. 

In  Section  III  we  discuss  the  results  of  an  application  of 
the  VFM  technique  to  a large  population  of  North  American 
earthquakes  and  explosions  at  the  Nevada  Test  Site  (NTS) 
recorded  at  the  19-element  short-period  Yellowknife  array 
in  Canada.  These  results  are  compared  to  previous  results 
obtained  for  Eurasian  events  recorded  at  LASA. 

In  Section  IV  is  a theoretical  analysis  which  addresses 
the  question:  under  what  circumstances  can  tectonic  release 

have  an  important  effect  on  the  teleseismic  short  period  P- 
wave  signature  of  underground  nuclear  explosions?  The 
analysis  has  validity  whether  one  assumes  that  the  primary 
mechanism  for  the  release  of  tectonic  stress  is  the  creation 
of  a spherical  shatter  zone  by  the  explosion  or  movement  along 
^ P^-e— exis ting  fault  plane.  In  either  case,  it  is  only  under 
optimal  conditions  that  the  tectonic  release  contribution  to 
the  seismogram  becomes  important.  It  is  concluded  that  in 
the  case  of  most,  if  not  all,  events  it  is  safe  to  ignore 
tectonic  release  as  a contributor  to  the  short  period  P-wave 
recording. 

A three-dimensional  finite  difference  simulation  of 
an  earthquake  is  discussed  in  Section  V.  The  calculation 
includes  a realistic  nonlinear  model  of  spontaneous  stick- 


slip  earthquake  faulting  for  the  case  of  a bilateral  rupture 
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equivalent  elastic  point  source  representation  of  the  earth 
quake  is  obtained.  Using  this  unique  and  exact  representa- 
tion, the  far-field  (that  which  propagates  to  teleseismic 
distances)  radiation  of  stress  waves  by  the  earthquake 
source  is  studied. 

Appendices  A - D pertain  to  key  elements  of  the  3-D 
finite  difference  earthquake  calculation  described  in  Sec- 
tion V.  Sections  II  and  III  should  be  read  in  sequence. 
Sections  IV  and  V are  self-contained  and  may  be  read  inde- 
pendently. 
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THE  MULTIPLE  ARRIVAL  RECOGNITION  SYSTEM 


INTRODUCTION 


Significant  progress  has  been  made  on  the  Multiple 
Arrival  Recognition  System  (MARS)  since  the  last  semi-annual 
report  (Bache,  et  al.  [1975a]),  where  it  was  briefly  outlined 
in  an  appendix.  As  will  be  shown  in  this  section,  the  MARS 
code  has  developed  into  an  extremely  powerful  and  efficient 
tool  for  signal  (body  and  surface  waves)  detection  and  en- 
hancement, and  for  discrimination  between  earthquaVes  and 
underground  explosions. 

improvements  made  to  the  MARS  code  during  the  last 
six  months  include  provisions  to  correct  envelopes  for  fre- 
quency-dependent noise,  to  remove  signal  distortion  due  to 
instrumental  factors,  and  to  perform  polarization  filtering 
on  multicomponent  data.  Improvements  have  also  been  imple- 
mented for  display  of  the  computed  results. 


NARROW-BAND  FILTERING  (NBF^. 


Data  are  input  to  MARS  in  the  form  of  a time  series, 
generally  of  about  500-2000  points  in  length.  The  data  are 
then  optionally  demeaned,  detrer.ded,  and  tapered  at  the  tail 
end  by  a cosine  bell.  MARS  then  ‘elects  the  smallest  power 
of  two  which  is  greater  than  the  number  of  points  input,  and 
performs  a discrete  Fourier  transform  using  the  algorithm  of 
Cooley  and  Tukey  [1965] . Both  the  time  series  and  the 
spectrum  are  plotted  for  examination. 

The  signal  is  filtered  in  the  frequency  domain  by 
multiplication  with  a cusp-type  filter  of  the  form: 
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if  - ( f_  - S-  Af  )' 


,f  - i Af  < f < f 
c 2 — — c 


F (f ) = < 1 - cos  — 


'(fc  + I A£)  - f 


Af 


, f < f < f + ^ Af 
c — - c 2 


0 , otherwise 


This  filter  is  drawn  in  Fig.  2.1.  It  was  discussed  in  the  pre- 
vious report  (Bache,  et  al.  [1975a] ) , compared  against  several 
other  filter  shapes,  and  found  to  be  the  best  in  terms  of 
time-domain  ripple  suppression.  The  width  at  one-half  maxi- 
mum amplitude  is  designated  Af. 

* I 


Frequency  ► 


Figure  2.1.  Filter  shape  used  in  MARS. 
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2.3  INSTRUMENTAL  CORRECTION 

Once  a signal  has  been  narrow-band  filtered,  MARS  cor- 
rects it  for  instrument  response.  When  a seismogram  is 
originally  taken,  it  is  measured  through  a system  whose 
transfer  function  is  strongly  frequency-dependent;  this 
distortion  affects  both  the  amplitude  and  the  phase  of  the 
signal.  If  we  are  to  make  accurate  measurements  of  magni- 
tude and  arrival  time,  the  signal  must  be  corrected  for  the 
appropriate  instrument  response. 

A typical  seismic  system  consists  of:  (a)  a trans- 

ducer, which  converts  ground  motion  into  electrical  signals; 
(b)  a galvanometer,  driving  a mirror  to  provide  optical 
amplification;  (c)  a system  of  photocells  to  provide  conver- 
sion back  to  an  analog  electrical  signal;  (d)  an  electrical 
filter  to  further  shape  the  pass-band,  and  finally,  (e)  a 
digitizer  and  recorder  where  the  data  are  stored. 

The  transducer  and  galvanometer  considered  separately, 
both  obey  second-order  ordinary  differential  equations.  When 
they  are  coupled  together,  the  resulting  system  obeys  the 
equation  (c.f.,  for  example,  Kisslinger  [1971]) 

( iv)  • • • v • ••• 

0 v ' + a e + b0  + ce  + dt  = -e  x 

for  the  galvanometer  deflection  angle  (0)  due  to  a ground 
motion  (X) , where 

A = ?hTu)T  + 2hGwG 

B - “ T + 4hThGuT“G  + “G  ■ 4hThG“TuG°2 

c " 21it“t“g  + 2Vg“t 

D = ^ 
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E = 


2hGU)G° 


coT  = transducer  resonant  frequency 


u)  = galvanometer  resonant  frequency 
G 


hT  = transducer  damping  fraction 


h_  = galvanometer  damping  fraction 
G 


a = a constant  which  measures  coupling  between 
transducer  and  galvanometer 


H = a measure  of  the  mass,  restoring  spring, 
and  level  arm  of  the  transducer 


Fourier  transformation  of  this  yields  the  transfer 
function 


X (“>)  = 


iw*/x. 


(u)4  - Bu>2  + D)  + i(-Aw3  + Cw) 


with  x being  chosen  as  a real  number  which  causes 
o 

|x  | = 1 at  a frequency  of  w/2tt  = f = 1 Hz.  An  approxima- 
tion to  a short-period  LRSM  instrument  is  given  by  using 
the  values 


= 2tt/1.02 


u)_  = 2tt/0 . 2 

b 


hT  = 0.98 


hG  = 0.9 


a2  = 0.01 


The  MARS  program  divides  the  signal  transform  by  this  transfer 
function,  after  the  former  has  been  narrow-ba  id  filtered. 
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The  electrical  filter  used  on  the  seismic  data  is 
typically  a Geotech  model  6824-1,  or  one  of  a similar  family. 
These  filters  can  be  described  reasonably  well  by  the  transfer 
function 

iwR2/X 

X(w)  ; . 

(a  f iw) [R2  - w2  + i2SRw] 


The  first  term  of  the  denominator  provides  a very  sharp  low- 
frequency  cut-off;  a is  generally  taken  to  be  tt/50.  The 
second  part  of  the  denominator  is  clearly  a harmonic  oscil- 
lator response;  for  this  we  have  taken  R = 40^/7  and 
S = /I72".  Like  the  instrument  response,  this  transfer 
function  is  normalized  to  be  1 at  1 Hz  by  a judicious  choice 
of  x . The  narrow-band  filtered  signal  transform  is  also 
divided  by  this  factor,  to  bring  the  signal  much  closer  to 
a filtered  true  ground  motion. 

The  resulting  complex  spectrum  is  then  inverse 
Fourier  transformed  into  the  time  domain,  to  produce  what 
will  hereinafter  be  called  the  filtered  signal. 

2.4  ENVELOPE  CONSTRUCTION  BY  HILBERT  TRANSFORMS 
(Bracewell,  [1965]) 

The  filtered  signal  produced  by  MARS  has  the  form 
X(t)  = A (t ) cos(wct  + 4>)  , wc  > 0 , 


which  clearly  represents  a modulated  carrier  wave.  Its 
Fourier  transform  is 


X(w) 


-/ 


X(t)  e"ia)t  dt  . 
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This  can  be  integrated  to  produce 

X(w)  = j [el4>  A(u>-u>c)  + e-1^  A<«+wc>]  , 

where  A(w)  is  the  Fourier  transform  of  A(t). 

If  we  consider  the  same  signal,  but  with  a 90°  phase 

shift, 

Y{t)  = A (t ) sin(oict  + <j>)  • 

then  its  transform  is 

Y(w)  = ~ [-iel4>  A(m-uc)  + ie_1<,>  A(w+«c)l  • 

These  two  functions  of  frequency  can  be  easily  related  by 
the  Hilbert  transform  condition, 

Y (w)  = -i  sgn(w)  X(w)  , 

provided  that  A(u>)  has  significant  amplitude  only  near 
w = 0,  and  drops  essentially  to  zero  for  M > uc-  This 
condition  is  equivalent  to  having  a narrow-band  signal; 
since  our  filtered  signal  obeys  this  condition,  we  can  gen- 
erate a signal  which  has  the  same  envelope  function,  A(t), 
as  the  filtered  signal,  but  which  is  90°  out  of  phase.  This 
secondary  signal  is  ordinarily  referred  to  as  the  quadra- 
ture, while  the  corresponding  analytic  signal  is  defined 

by 

i(o3ct  + 4)) 

Z(t)  = X(t)  + i V(t)  = A (t ) e 

It  is  of  interest  to  note  that  the  analytic  signal  contains 
only  positive  frequency  components.  The  envelope  is  ex- 
tracted from  it  quite  easily  as 

A(t)  = | Z It)  | = Vx‘  (t)  + Y2  (t)  , 
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while  the  instantaneous  phf.se  and  frequency  are  given  as 

<t(t)  = w t + if>  = Arctan[Y(t)/X(t)  ] , 

c 

and 


u 


O 


o 


Q 


U)(t)  = 


d$ 

dt 


= U) 


d$ 

dt 


This  method  is  followed  in  MARS:  The  transform  of  the 

filtered  signal  is  multiplied  by  -i  sgn (w) , and  then  brought 
to  the  time  domain  by  an  inverse  transformation.  The  envelope 
and  instantaneous  frequency  are  produced;  the  envelope  has 
applications  in  terms  of  variable-frequency  magnitude  (VFM) 
measurements,  while  the  instantaneous  frequency  is  utilized 
mostly  for  polarization  filtering.  Both  of  these  subjects 
will  be  discussed  later  and  at  greater  length. 

An  envelope,  constructed  by  narrow-band  filtering  in 
this  manner,  will  be  a very  smooth  function  of  time.  The 
time  required  for  the  envelope  to  change  its  height  appreci- 
ably must  be  proportional  to  the  inverse  filter  width  1/Af  = 
Q/f ; this  statement  is  essentially  just  the  well-known  un- 
certainty principle.  Stated  in  terms  of  information  theory, 
a signal  of  a particular  bandwidth  is  limited  in  the  amount 
of  information  which  it  can  carry,  so  that  the  highest  rate 
of  information  flow  is  the  same  as  the  frequency  width  of 
the  channel . 

it  is  true  that  a signal  of  a pure  frequency  has  no 
time  dependence.  However,  it  is  possible,  using  only  those 
frequency  components  which  fall  within  a limited  bandpass, 
to  make  up  a wave-packet  group  which  is  fairly  localized  in 
time.  Taken  in  this  light,  the  time  variation  of  an  enve- 
lope can  be  interpreted  as  being  the  arrival  (or  non-arrival) 
of  wave  groups  of  the  filter  frequency;  the  variation  will 
describe  the  appearance  of  energy  of  that  frequency. 
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2.5 


DETECTION  OF  ARRIVALS  IN  THE  t 

a 


f PLANE 
c 


Given  a particular  seismic  signal,  the  procedure  of 
envelope  construction  can  be  performed  for  a number  of  fil- 
ter center  frequencies,  fc«  For  each  of  these  envelopes,  a 
group  arrival  time  (tg)  can  be  computed  and  plotted  as  a func' 
ti.on  of  f t.  If  the  original  signal  consists  of  only  one  un- 
dispersed  arrival,  then  the  arrival  time  diagram  will  look 
like  Fig.  2.2a. 

If  a single  arrival  is  present,  but  it  has  been  dis- 
persed such  that  low  frequencies  travel  fastest  (e.g., 
normally  dispersed  Rayleigh  waves)  then  the  arrival  time  dia- 
gram looks  like  Fig.  2.2b.  Figure  2.2c  exhibits  the  opposite 
effect,  with  high  frequencies  traveling  faster  than  lows,  as 
in  the  case  of  inversely  dispersed  surface  waves.  Under  cer- 
tain circumstances,  a dispersion  curve  can  bend  down  and  have 
a flat  portion,  oc  a minimum.  This  type  of  dispersion  gives 
rise  to  an  Airy  phase  (B&th  [1968])  and  is  demonstrated  in 
Fig.  2. 2d. 


2.6  VARIABLE  FREQUENCY  MAGNITUDE  DETERMINATION 


For  each  filter  frequency  and  for  each  significant 
signal  arrival,  in  terms  of  signal-to-ncise  ratio,  it  is  pos- 
sible to  compute  a seismic  magnitude.  The  procedure  followed 
in  MARS  is  to  define  body  wave  magnitudes  as 


m.  = log  (S/f  ) + B (A)  + 0.05 
b io  max  c 

and  surface  wave  magnitudes  as 

Ms  - lo910(Wfc>  + i-656  l°910(4>  + l-6  • 
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Figure  2.2a.  Arrival  time 
diagram  for  an  undispersed 
pulse. 


Figure  2.2b.  Arrival  time 
diagram  for  a pulse  which  has 
been  dispersed  by  propagation 
through  a medium  in  which  low 
frequencies  travel  faster  than 
high. 


Figure  2.2c.  Dispersion  causes 
high  frequencies  to  travel 
fastest. 


Figure  2. 2d.  Airy  phase  dis- 
persion. 
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distance  between  event  and  detector  in 
decrees. 


The  function  B(£)  represents  modified  Gutenberg-Richter 
distance  correction  factors,  derived  from  LRSM  data.  The 
exact  values  used  are  shown  in  Table  2.1. 


TABLE  2.1.  MODIFIED  GUTENBERG-RICHTER  B-FACTORS 
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2.7  NOISE  CORRECTION 

When  an  event  is  processed  which  contains  a signifi- 
cant amount  of  noise,  a noise  correction  can  be  made  to  the 
spectral  amplitudes.  This  is  done  by  taking  a section  of 
the  seismogram  before  the  arrival  of  the  actual  event,  and 
processing  the  noise  alone  as  though  it  were  a signal.  The 
noise  envelope  is  formed  for  each  of  several  filter  frequen- 
cies, and  in  each  case  the  maximum  value  is  taken  as  an 
amplitude  correction  factor. 

Subsequently,  when  the  signal  window  is  processed  and 
the  signal  envelopes  constructed,  the  maximum  noise  amplitude 
at  that  frequency  is  subtracted  as  a constant.  This  will 
cause  the  less  significant  signal  envelope  peaks  to  be  sub- 
merged, and  they  can  be  disregarded  as  far  as  arrival  identi- 
fication is  concerned. 

2.8  DEMONSTRATION  OF  NARROW-BAND  FILTERING  FOR  SEISMIC 

EVENTS 

In  this  subsection  of  the  report  the  narrow— band 
filtering  technique  described  above  will  be  applied  to  actual 
seismic  signals  recorded  at  the  LASA  and  NORSAR  short-period 
vertical-component  arrays  located  in  Montana  and  Norway, 
respectively.  The  seismograms  are  best  beam  recordings  based 
on  phasing  of  the  full  arrays. 

Figure  2.3a  shows  a LASA  best  beam  recording  of  a 
P-wave  train  from  a presumed  explosion  that  occurred  in 
eastern  Kazakhstan:  hereafter  referred  to  as  LASA  Event  1522 

The  amplitude  scale  given  along  the  ordinate  in  this  figure 
is  in  millimicrons  of  grourd  motion  at  1 Hz.  Figure  2.3b  is 
the  amplitude  spectrum  of  :his  event,  showing  that  the  signal 
energy  is  principally  in  the  region  of  1 Hz,  as  one  would 
expect  for  body  waves  detected  at  teleseismic  distances  by 
short-period  instruments  of  the  LASA  type. 
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The  time  series  in  Fig.  2.3a  was  processed  through  ten 
narrow-band  filters,  with  center  frequencies  tQ  = 0.4,  0.6, 

0.8,  1.0,  1.2,  1.4,  1.5,  1.6,  1.8  and  2.0  IIz,  all  with  quality 
factors  of  Q = 10.  The  single  arrival  nature  of  these  data 
is  reflected  in  Fig.  2.3c,  which  shows  the  sun  of  the  envelopes 
for  the  ten  filters  as  a function  of  tine. 


Figure  2.3d  is  an  arrival  time  diagram,  on  which  time 
runs  horizontally  and  filter  center  frequency  vertically.  For 
each  frequency,  the  envelope  is  plotted  as  a function  of  time; 
each  of  these  curves  has  been  normalized  so  as  to  have  the 
same  height.  The  undispersed  nature  of  these  data  are  quite 


evident. 


A LASA  recording  of  a P-wave  train  from  a Japanese 
earthquake  is  shown  in  Fig.  2.4a;  this  is  referred  to  as 
LASA  Event  2030.  In  it,  a direct  P-wave  is  demonstrated 
(arriving  at  approximately  t = 20  seconds)  followed  by  a sur- 
face reflection  (pP  at  * t = 31  seconds) . This  event  was 
processed  through  the  same  filters  as  for  the  previous  event 
(LASA  1522),  and  the  sum  of  the  envelopes  is  shown  in  Fig. 

2.4b.  The  separation  between  the  envelope  peaks  is  11.0  sec- 
onds, compared  with  the  separation  of  approximately  12.0  sec- 
onds in  the  unfiltered  signal.  This  event  was  reported  in  the 
Preliminary  Determination  of  Epicenters  (PDE)  to  be  at  a depth 
of  40  km  and  a distance  of  63.5°  from  LASA,  for  which  a pP-P 
time  delay  of  11.4  seconds  is  considered  normal  (Herrin, 
et  al . [1968] ) . 

For  Event  2030,  the  envelopes  for  different  frequen- 
cies are  shown  in  Fig.  2.4c.  This  figure  demonstrates  an 
interesting  phenomenon:  at  low  frequencies,  an  envelope  maxi- 

mum appears  in  the  neighborhood  of  t = 32-33  seconds,  but 
vanishing  when  fQ  > 1 Hz.  At  the  same  time,  the  envelope  peak 
near  t = 22-23  seconds  grows  in  significance  with  increasing 
frequency.  The  message  conveyed  by  this  is  clea''  - the  plJ 
signal  does  not  contain  as  much  high  frequency  energy  as  the 
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Figure  2.4a.  LASA  Event  2030:  Figure  2.4b.  LASA  Event  2030 

Vertical  component.  Sum  of  ten  envelopes. 


Envelopes  for  Z -Component 


0.0  I 1 
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Figure  2.4c.  Envelopes  for  LASA  2030  vertical 
motion. 
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P signal.  This  can  possibly  be  interpreted  as  the  effect  of 
absorption  of  the  high  frequency  components  of  the  free  sur- 
face reflection,  pP,  by  relatively  low  Q,  high  attenuation, 
surface  l^'ers  in  the  source  region. 

bASA  Event  2011,  a Rumanian  earthquake,  is  shown  in 
Fig.  2.5a.  This  signal  was  filtered  at  the  same  frequencies 
as  the  previous  two  events,  with  an  arrival  time  diagram  of 
the  sum  of  the  envelopes  being  shown  in  Fig.  2.5b  and  one  for 
the  individual  filter  envelopes  in  Fig.  2.5c.  The  largest 
peak,  at  t = 23.0  seconds,  is  due  to  direct  P-wave;  the  next 
largest  is  at  t = 62.0  seconds.  This  delay  of  39.0  seconds, 
for  a source  at  a distance  of  76.2°  and  a depth  of  158  km 
(PDE) , falls  in  the  range  given  for  pP  in  the  travel-time 
tables  of  Herrin,  et  al.  [1968]  whose  delay  times  are  repro- 
duced in  Table  2.2. 

TABLE  2.2.  pP  DELAY  TIMES  IN  SECONDS 


Distance 


Depth 


150  km  200  km 


36.2 


47.0 


36.6 


47.6 


Examination  of  Fig.  2.5c  shows  that  the  amplitude  of  the  pP 
phase  dies  out  with  increasing  filter  frequency,  relative  to 
the  amplitude  of  the  direct  P-wave:  below  1 Hz,  the  pP  is 

as  strong  or  stronger  than  the  P,  while  above  that  frequency, 
direct  P is  clearly  dominant.  This  is  similar  to  the  be- 
havior as  noted  for  Event  2030. 

The  third  largest  peak  in  Fig.  2.5b  lies  at  t = 33.0 
seconds,  a delay  with  respect  to  P of  10.0  seconds.  This 
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fits  in  well  with  Herrin's  values  for  a core  reflection,  PcP, 
and  the  bracketing  data  are  depicted  in  Table  2.3. 

TABLE  2.3.  DELAY  TIMES  BETWEEN  P AND  PcP 


Distance 


Depth 


150  km  200  km 


13.0 


12.7 


Figure  2.5c  shows,  interestingly  enough,  that  the  PcP  signal 
strength  is  roughly  a constant  fraction  of  the  direct  P 
amplitude.  Note  that  as  we  filter  at  higher  and  higher  fre- 
quencies, the  PcP  phase  remains  visible  with  about  20-30  per- 
cent of  the  maximum  envelope  height  associated  wit  the  P 
phase.  This  shows  that  PcP  has  not  suffered  the  high-frequency 

attenuation  observed  for  pP. 

On  Fig.  2.5b,  a fourth  peak  is  visible  near  t = 42-43 
seconds.  Examination  of  Fig.  2.5c  indicates,  however,  that 
the  frequency  content  of  this  particular  peak  is  within  the 
microseisrdc  band. 

The  seismogram  for  a P-wavetrain  from  an  earthquake  in 
the  Caucasus  recorded  at  NORSAR  (NORSAR  Event  6277),  is  shown 
in  Fig.  2.6a.  It  is  clearly  a very  noisy  record.  This 
seismogram  was  processed  by  MARS,  using  the  previously 
described  method  designed  to  eliminate  noise. 

The  noise  window  was  taken  to  be  the  section  of  the 
recording  in  Fig.  2.6a  from  t = 0 to  t = 20  seconds,  the 
spectrum  of  this  portion  is  shown  in  Fig.  2.6b.  The  signal 
window  was  defined  as  the  51.2  second  interval  (1024  points 
at  20  points/second)  following  t = 20  seconds.  The  signal 
spectrum  is  plotted  in  Fig.  2.6c.  Comparison  of  these  two 
spectra  makes  clear  the  profound  difference  in  frequency 
content  between  noise  and  signa.'. . 
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When  the  sum  of  envelopes  is  formed,  with  the  appro- 
priate noise  correction  subtracted  from  each  individual  enve- 
lope, the  result  is  as  shown  in  Fig.  2.6d.  In  this  figure 
the  straight  line  at  zero  amplitude  indicates  the  noise  level 
and  eliminates  all  but  two  peaks.  In  this  plot,  the  time 
scale  has  been  shifted  by  an  amount  such  that  the  signal 
window  begins  at  t = 40  seconds.  The  largest  peak,  at  t = 

43  seconds,  is  then  the  direct  P-wave  arrival;  the  second 
peak  is  about  9 seconds  later  - a plausible  delay  time  for 
pP,  if  the  event  was  in  the  neighborhood  of  35  km  deep. 

The  individual  envelopes  for  Event  6277  are  shown  in 
Fig.  2.6e,  with  the  noise  subtracted,  an  operation  which 
causes  the  lowest  frequency  envelope  (0.4  Hz)  to  virtually 
disappear  due  to  the  large  amount  of  long-period  background. 
In  this  display,  the  second  arrival  is  not  definitely  indi- 
cated (it  can  best  be  seen  in  the  f = 1.2  trace),  showing 
that  both  types  of  display  are  necessary  for  a good  signal 
evaluation. 

2.9  POLARIZATION  FILTERING 

Given  seismograms  with  more  than  one  component  of 
motion,  it  is  possible  to  perform  additional  analysis  based 
upon  the  phase  correlation  of  the  filtered  particle  motion. 

It  is  well-known,  for  example,  that  if  e is  a unit  vector 
oriented  vertically,  and  er  is  a unit  vector  pointing  away 
from  the  source,  the  P-wave  motion  will  tend  to  appear  as 

X(t)  = f(t)  [e„  cosa  + e„  sina]  , 

•«*  r z 

where  a is  the  apparent  emergence  angle.  SV  waves  will  be 
given  by 

X (t)  = g(t)  (-ev  cosa  + e„  sina] 

i Z 
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Figure  2.6e. 


Individual  envelopes  for  Norsar  6277 
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Both  of  these  are  linearly  polarized  in  the  z-r  plane,  the 
first  being  in  phase,  and  the  second,  180°  out  of  phase. 

Rayleigh  waves,  on  the  other  hand,  are  elliptically 
polarized,  having  the  form 

X ( t ) = e a cos<J>  (t)  + e b s±n4>  (t)  , 

«.  r « 

so  that  the  two  components  are  90°  out  of  phase  with  one 
another.  If  the  amplitude  coefficients  a and  b are 
equal,  then  the  polarization  is  circular,  and  Rayleigh 
wave  motion  is  often  close  to  this  limit. 

If  both  z-  and  r-component  motion  are  measured  simul- 
taneously, narrow-band  filtering  can  be  performed  on  the 
two  for  the  same  filter  frequencies  and  then,  for  each  ..'re- 
quency,  the  linear  polarization  fraction  can  be  defined  using 
the  instantaneous  phases  as 


PL(t)  = cos[(j)z{t)  - 4>r(t)]  , 

and  the  circular  polarization  fraction  as 
Pc(t)  - sin[<J>z  (t)  - 4»r  (t)  ] 


Clearly, 


K + pc  = 1 


so  that  the  two  are,  in  a sense,  mutually  exclusive. 

This  method  can  be  used  to  identify  different  signals 
as  belonging  to  various  wave  types.  Since  frequency  filter- 
ing is  done  before  polarization  filtering,  the  normal  sort  of 
low-frequency  Rayleig’,  wsve  background  noise  can  be  excluded 
at  an  early  stage,  allowing  the  detection  of  higher  frequency 


O linearly-polarized  information  with  greater  sensitivity.  If, 

I 


on  the  other  hand,  a P-vave  and  an  SV-wave  arrive  . imultan- 
eously,  then  they  may  tie  indistinguishable  from  a higher-mode 


Rayleigh  wave. 

For  z-r  data,  different  types  of  waves  can  be  identi- 
fied as  follows: 

p-wave : PL  > 0 and  P^  > I I 

SV-wave:  PL  < 0 and  |PL|  > I I 

Rayleigh  wave:  Pc  < 0 and  jPc|  > I PL l 

The  case  of  Pc  > 0 and  | PQ | > |PL!  can  be  taken  to  indicate 
either  noise,  or  higher  mode  Rayleigh  motion. 

For  waves  vhich  appear  to  be  linearly  polarized,  an 
apparent  emergence  angle  can  be  defined  by 

a(t)  = Arctan[Sz (t)/Sr (t) ] , 

while  the  ellipticity  of  other  waves  can  be  defined  by 

E (t)  = Sz(t)/Sr(t)  . 

If  the  data  have  not  been  rotated  exactly  towards  the 
source,  and  if  three  Cartesian  components  are  available, 
then  the  same  analysis  can  be  performed,  comparing  the  verti- 
cal with  each  of  the  horizontal  components  separately.  Cor- 
relation of  the  two  analyses  can  then  be  used  to  indicate  the 
apparent  direction  of  the  source;  since  differing  frequency 
components  will  often  be  refracted  differently,  this  apparent 
direction  may  be  expected  to  vary  somewhat  with  frequency. 

Three-component  data  can  also  be  used  to  identifv 
Love  waves.  If  the  data  have  already  been  analyticall 
rotated  toward  the  source,  then  the  transverse  channel  will 
contain  mostly  Love  and  SH  waves.  Should  the  signal  not 
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have  been  rotated,  then  these  wave  types  will  be  linearly 
polarized  in  the  horizontal  plane,  with  an  angle  which  will 

indicate  an  apparent  azimuth. 


2.10  CONSTRUCTION  OF  POLARIZATION-FILTERED  WAVES 


A two-dimensional  space  can  be  constructed,  the  basis 
set  of  which  consists  of  one  linear  and  one  circular  component. 
The  linear  and  circular  polarization  fractions  for  a given 
wave  then  indicate  where  the  signal  appears  in  such  a*  space. 
The  polarization  fractions  can  be  used  to  form  projections 
into  this  space,  where  individual  wave  types  can  be  easily 
picked  out.  For  example,  the  P-wave  envelope  projection  can 
be  defined  as 


Sp  (t)  - max(PL(t)  ,0)  • VsJ(t)  + Sj(t)  , 


the  SV-wave  projection  as 

Ssv(t)  = max(-PL(t)  ,0)  * Vs;  (t)  + S*  (t)  , 


and  the  Rayleigh  wave  projection  as 


SR(t)  - max(-Pc(t),0)  • VsJ(t)  + Sj(t)  . 


By  evaluating  these  quantities  for  several  frequencies,  the 
appearance  of  different  wave  types  can  be  observed  and  separated 

by  polarization. 


2.11  DEMONSTRATION  OF  POLARIZATION  FILTERING  FOR  ARTIFICIAL 
SIGNALS 


To  test  the  polarization  filtering  method,  an  artifi- 
cial two— component  signal  was  prepared,  with  the  intention 
that  it  look  like  a P-wave,  followed  by  an  SV-wave.  Both 
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were  represented  by  decaying  sine-waves,  the  P-wave  as 
z(t)  = r (t)  = 0<t-15)e"(t"15)/5  sin  (t-15) 


and  the  SV  as 


z(t)  -=  -r  (t)  = 0.4  0 (t-25)e“(t“25)/4  sin 


t-25 


where  0 is  the  Heaviside  step  function.  The  two  waves 
can  be  seen  to  have  slightly  different  frequency  content,  as 
well  as  different  decay  times:  The  P-wave  train  should  last 

a bit  longer  than  the  SV. 

The  z -component  is  shown  in  Fig.  2.7a,  and  the  r in 
Fig.  2.7b.  The  SV-wave  onset  can  be  seen  in  the  first  dia- 
gram by  a sharp  break  at  t = 25;  it  is  not  visible  in  the 
second  plot,  and  it  is  not  obvious  that  this  is  anything  but 
a single  complicated  wave.  About  1 1/2  cycles  separate  the 
two  waves,  and  so  the  P— wave  will  have  to  be  observed  during 
this  time.  In  addition,  detection  of  the  SV-wave  will  be 
hampered  by  the  coda  of  the  first  arrival,  which  has  a larger 
initial  amplitude  and  greater  decay  time. 

This  time-series  (P  and  SV)  was  filtered  by  ten  filters 
with  center  frequencies  between  0.25  and  3.0.  When  the  enve- 
lopes are  summed,  the  result  is  as  shown  in  Figs.  2.7c  (for  z) 
and  2.7d  (for  r) . In  both  these  cases,  the  SV  arrival  stands 
out  clearly.  The  peaks  for  both  arrivals,  on  both  components, 
match  the  actual  onset  times  qul^e  accurately. 

The  individual  envelopes  are  plotted  in  Figs.  2.7e 
and  2 . 7f . Here  we  see  that  the  narrow-band  filtering  has 
done  a very  good  job  of  separating  the  arrivals  at  the  higher 
frequencies. 

When  the  circular  polarization  fraction  (Pc)  is  cal- 
culated for  each  filter  center  frequency,  as  a function  of 
time,  Fig.  2.7g,  the  results  are  essentially  nil.  Aside  from 


32 


Figure  2.7a.  Vertical  component  of  Figure  2.7b.  Radial  component  of 

displacement.  displacement. 
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a spike  near  t = 20  (where  it  separates  the  P and  SV  waves) , 
no  significant  circular  arrivals  can  be  seen. 

The  linear  polarization  fraction  (PL)  , in  Fig.  2.7h, 
on  the  other  hand,  indicates  the  two  wave  types  with  great 
clarity.  It  can  be  seen  that  there  is  a nearly  constant 
positive  value  (in  fact,  +1)  at  early  times,  switching 
rapidly  to  a value  below  zero  (generally,  -1)  at  t = 25. 

At  much  later  times  (t  = 40-60),  the  lower  damping  of  the 
P-wave  makes  itself  felt,  and  its  coda  begins  to  dominate 
that  of  the  shear  excitation. 

■"lie  Rayleigh  wave  projections  (SR)  are  plotted  in 
Fig.  2.7i  as  functions  of  time,  for  the  various  filter  fre- 
quencies. These  data  are  quite  inconsistent,  every  other 
filter  indicating  a different  motion. 

The  P-wave  (Sp,  Fig.  2 . 7 j ) and  SV-wave  projections 
(Sgv,  Fig.  2.7k)  show  very  obvious  arrivals.  The  lower 
resonant  frequency  of  the  P-wave  results  in  the  apparent  de- 
lay of  the  SV-wave  envelope  at  the  lowest  frequency  - much 
more  P-wave  energy  than  SV  arrives  for  this  frequency,  so 
that  the  latter  has  trouble  making  its  appearance. 

This  test  was  then  rerun,  with  a shorter  delay  time  be- 
tween the  arrivals  — the  SV  wave  arrival  time  was  advanced  from 
t = 25  to  t = 20  seconds.  The  z-component  is  shown  in  Fig.  2.8a, 
and  the  r in  Fig.  2.8b.  The  second  arrival  is  not  easily 
seen  in  either  of  these  records;  the  narrow-band  filtering  sum- 
of-envelopes  plots  (Figs.  2.8c  and  2.8d)  reveal  its  presence  by 
a slight  bump.  Individual  envelopes,  in  Figs.  2.8e  (z-component) 
and  2.8f  (r-component ) show  the  shear  wave  at  high  filter  fre- 
quencies. By  looking  at  Fig.  2.8g,  we  can  see  the  sane  strange 
behavior  of  the  circular  fraction  as  was  observed  earlier;  the 
linear  fraction,  in  Fig.  2.8h,  however,  still  shows  the  very 
strong  identification  of  the  two  waves  with  their  different 
polarizations.  The  projected  P-wave  envelopes  of  Fig.  2.3i, 
and  SV-wave  envelopes  of  Fig.  2 . 8 j are  almost  as  they  wore  in 
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Figure  ?.7i.  Rayleigh  wave  pro- 
jection, based  upon  circular  polar- 
ization fraction. 


Figure  2.7j.  P-wave  projection 
as  obtained  from  linear  polar- 
ization fraction. 
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Figure  2.7k.  SV-wave  projection,  as  ob- 
tained from  linear  polarization  fraction. 
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the  first  case,  demonstrating  that  the  polarization  analysis 
is  still  successful  in  locating  the  arrivals. 


The  test  was  then  toughened  by  moving  the  shear  wave 
forward  some  more:  The  P-wave  now  begins  at  t = 15,  and  the 

SV  at  t = 17.5.  Figures  2.9a  and  2.9b  display  the  z and 
r displacements,  respectively.  No  indication  of  the 
multiple-arrival  nature  can  be  seen  by  the  eye  on  either  of 
these  plots.  The  narrow  band  discriminant  fails,  as  the  sum 
of  z-envelopes,  Fig.  2.9c,  shows;  the  individual  z-envelopes, 
Fig.  2.9d,  indicate  the  second  arrival  clearly  only  at  the 
highest  frequency,  although  some  of  the  lower  frequencies 
exhibit  a prolonged  decline  where  the  SV-wave  should  be. 


The  linear  polarization  fraction  (Fig.  2.9e)  is  shaken, 
but  not  fooled.  The  SV-arrival  is  indicated  plainly  on  three 
of  the  traces  (f  = 1.0,  2.0,  and  3.0)  and  is  hinted  at  by 
the  other  curves. 


The  P-waves  are  projected  and  drawn  in  Fig.  2.9f;  their 
presence  is  evident.  The  appearance  of  SV-waves  is  also  easily 
demonstrated  in  Fig.  2.9g,  although  again  the  lower  frequencies 
suffer  from  P-wave  domination. 


2.12  POLARIZATION  ANALYSIS  OF  SEVERAL  OBSERVATIONAL  EVENTS 


Figures  2.10a  and  2.10b  show  vertical,  z,  and  radial, 
r,  component  seismograms,  respectively,  of  ground  motion  from 
the  Rex  explosion  at  the  Nevada  Test  Site  (NTS)  recorded  at 
the  LRSM  station,  RKON.  These  two  records  were  processed  by 
MARS  for  filters  between  1 and  4 Hz,  resulting  in  the  enve- 
lopes shown  in  i^igs.  2.10c  and  2.10d.  The  displacement  in 


the  z-r  plane  waj  decomposed  into  circular  (P  ) and  linear  (PT ) 

C Ij 


fractions,  as  exhibited  in  Figs.  2.10e  and  2.10f. 


A comparison  between  Figs.  2.10c,  2.10e,  and  2.10f  reveal 
this  fact:  When  the  envelopes  are  large,  the  circular  fraction 

tends  to  be  near  zero,  while  the  linear  fraction  is  most  likely 
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Figure  2.9e.  Linear  polarization 
fraction  for  the  Z-R  plane.  Note 
that  a change  in  polarization  is 
consistently  indicated  near 
t = 17.5. 
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Figure  2.9f.  Projected  P-waves 
using  linear  polarization  frac- 
tion of  Fig.  2.9e. 


foUili<tiOA-mtir<4  biultn  SV-ttavaa 


0.0  20.0  40.0  (0.0  10.0  100.0  120.0 

TlM  doc) 

Figure  2.9g.  SV-waves  projection,  using  linear 
polarization  fraction  of  Fig.  2.9e.  Note  that, 
for  higher  frequencies,  the  onset  time  is  con- 
sistent. 
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to  be  near  +1  and  flat.  This  is  a very  clear  indication  that 
the  observed  excitation  is  predominantly  of  a longitudinal 
(P-wave)  nature. 

When  we  attempt  to  construct  the  Rayleigh  wave  envelope 
projection  (SR) , the  results  are  ragged  and  inconsistent,  as 
shown  in  Fig.  2.10g.  A comparison  of  these  traces  with  the 
envelope  amplitudes  show  that  Rayleigh  waves,  likely  compris- 
ing the  microseismic  background,  are  indicated  generally  at 
times  when  no  signal  energy  i«*  arriving.  It  should  be  noted, 
in  viewing  this  diagram,  that  each  of  the  curves  has  been 
magnified  so  that  its  maximum  value  shows  as  a standard  size 
on  the  plot;  the  relative  lack  of  strength  of  this  wave  type 
is  therefore  not  demonstrated. 

Figure  2.10h  shows  the  arrival  of  P-waves,  as  calculated 
by  the  P-wave  projection  (Sp) . A comparison  of  this  with  Fig. 
2.10c  reveals  that  most  of  the  envelope  pv.lses  are  of  this 
species,  a fact  which  was  noted  earl.  by  observation  of  the 
polarization  fractions.  Figure  2>10i  demonstrates  the  SV-waves, 
and  this  shows  that  the  record  is  pretty  clearly  devoid  of  such 
excitations. 

Figure  2.11a  is  the  z-component  of  the  NTS  explosion, 
Greely,  as  measured  at  RKON . The  corresponding  r-component  is 
shown  in  Fig.  2.11b  while  the  z and  r envelopes  are  in 
Figs.  2.11c  and  2. lid,  respectively.  The  circular  polariza- 
tion fraction  (Fig.  2. lie)  is  basically  uninteresting:  When 

the  envelopes  show  arrivals,  this  quantity  is  near  zero,  and 
nowhere  can  we  see  the  broad,  flat  areas  which  indicate  the 
arrival  of  a coherent  wave.  The  linear  polarization  fraction 
(Fig.  2. Ilf),  on  the  other  hand,  is  just  what  one  should  expect 
if  a definite  wave  type  has  arrived.  The  linear  part  is  almost 
+1  at  all  times  of  significant  arrival;  the  conclusion  is  that, 
without  a doubt,  we  are  looking  at  a P-wave  arrival. 
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Figure  2.10g.  Rex  at  RKON 
Rayleigh  wave  projections. 


Figure  2.10i.  Rex  at  RKON 
SV-wave  projections. 
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Figure  2.10h.  ReX  at  RKON:  P-wave  projections 
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The  projected  Rayleigh  (Fig.  2.11g) , SV  (Fig.  2.11h) , 
and  P waves  (Fig.  2. Hi)  verify  that  result.  All  of  the  sig- 
nificant peaks  in  the  z-envelope  have  been  projected  as  longi- 
tudinal rectilinear  waves. 

The  Piledriver  explosion,  recorded  on  the  long-period 
instruments  at  RKON,  is  shown  in  Fig.  2.12a  (z-component)  and 
2.12b  (r-component) . The  envelopes  due  to  filtering  at  dif- 
ferent frequencies  are  in  Figs.  2.12c  and  2.12a.  The  arrival 
time  is  manifestly  different  for  different  center  frequencies, 
with  the  low  frequency  energy  arriving  first.  The  variation 
of  peak  arrival  time  with  frequency  is  smooth  and  continuous, 
making  it  evident  that  this  is  due  to  a single  arrival,  rather 
than  several  overlapping  arrivals.  This  is  obviously  an  Airy 
phase. 

The  circular  polarization  fraction  (Fig.  2.12e)  is 
also  easy  to  interpret:  The  envelope  peaks  correspond  to 

areas  of  strong  and  consistent  negative  circularity  — this  is 
a Rayleigh  wave.  The  linear  fraction,  Fig.  2.12f,  on  the  other 
hand,  avoids  the  peaks.  In  Fig.  2.12g,  the  Rayleigh  waves  have 
been  projected  out,  and  the  Airy  phase  dispersal  is  again 
plainly  demonstrated. 

Figure  2.13a  is  the  transverse  component  of  the  same 
event.  A Love  wave  can  be  observed,  which  arrives  earlier 
than  the  Rayleigh.  The  arrival  time  diagram  for  this  (Fig. 
2.13b)  also  shows  the  unusual  pattern  characteristic  of  an 
Airy  effect. 
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Figure  2.11g.  Greely  at  RKON : Figure  2.11h.  Greely  at  RKON 

Rayleigh  wave  projections.  SV-wave  projections. 
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Figure  2. Hi.  Greely  at  RKON:  P-wave  projections. 
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III.  VARIABLE  FREQUENCY  MAGNITUDE  DISCRIMINANT 
3.1  INTRODUCTION 

In  a previous  report  (Bache,  et  al . , 197 Sa)  written 
under  this  contract,  the  formulation  and  results  of  exten- 
sive testing  of  a variable  frequency  magnitude  (VFM)  dis- 
criminant which  exploits  spectral  differences  between  earth- 
quakes and  underground  explosions  were  presented  in  detail. 
In  addition  to  providing  discrimination  of  a large  popula- 
tion of  shallow  Eurasian  earthquakes  and  presumed  explosions 
recorded  at  the  LASA  and  NORSAR  arrays,  it  was  shown  that 
the  VFM  technique  can  discriminate  a multiple  explosion 
scenario  as  well.  A particularly  significant  aspect  of  the 
VFM  technique  is  that  its  application  is  dependent  on  re 
cording  short-period  P-waves  from  events  and  not  long- 
period  surface  waves.  Thus,  positive  discrimination  can  be 
tested  for  on  events  (explosions)  down  to  a much  lower 
magnitude  threshold  level  than  is  possible  with  the  tradi- 
tional Mg  - m^  method. 

In  the  following  sub-sections  of  this  report  we  will 
present  preliminary  results  on  the  application  of  the  VFM 
technique  to  a large  population  of  North  American  events 
recorded  in  Canada.  This  section  will  conclude  with  a 
comparison  of  all  the  VFM  discrimination  tests  conducted 
to  date. 


3 . 2 NORTH  AMERICAN  EVENTS  RECORDED  IN  CANADA 

Data  in  digital  format  for  a large  population  of  North 
American  earthquakes  and  underground  explosions  recorded  at 
the  19-element  short-period  Yellowknife  array  in  Canada 
(Manchee  and  Somers,  1966)  were  recently  acquired.  These 
data  consist  of  vertical-component  recordings  of  P-wave 
trains  (preceded  by  at  least  25  seconds  of  background  noise) 
from  approximately  40  earthquakes  occurring  in  the  Gulf  of 
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California  and  49  underground  explosions  detonated  at  the 
Nevada  Test  Site  (NTS) . While  signals  recorded  on  all  19 
array  elements  were  supplied  for  most  of  the  events  in  the 
data  base,  our  first  test  for  discrimination  with  the  VFM 
technique  was  performed  on  single-element  (sensor  B5  located 
at  the  crossover  point  of  the  L-shaped  array  arms)  record- 
ings only. 

Examples  of  P-wave  trains  from  two  NTS  explosions  and 
two  Gulf  of  California  earthquakes  recorded  by  the  single 
sensor  B5  of  the  19-element  Yellowknife  array  are  shown  in 
Figs.  3. la-3. Id.  These  particular  time  series  were  chosen 
for  display  to  give  some  indication  of  the  variation  in 
signal- to-noise  ratio  (S/N)  of  the  different  event  signals 
in  the  data  base.  For  instance,  the  S/N  for  the  NTS  explo- 
sion Wagtail  (Fig.  3.1a)  is  more  than  a factor  of  ten  higher 
than  the  S/N  for  the  explosion  Diana  Mist  (Fig.  3.1b). 

Each  of  the  time  series  tested  for  discrimination  in 
this  study  was  first  divided  into  a noise  sample  and  a 
signal  sample.  In  the  case  of  the  explosions,  both  noise 
and  signal  windows  were  chosen  to  be  24  seconds  long.  For 
the  earthquakes,  a 24  second  noise  window  and  signal  win- 
dows up  to  66  seconds  in  duration  were  examined.  The  longer 
duration  of  the  earthquake  signal  windows  was  dictated  be- 
cause of  the  difficulty  in  identifying  the  exact  signal 
onset  time.  Each  noise  and  signal  window  was  demeaned  and 
tapered  with  a cosine  bell  applied  to  10  percent  of  both 
window  ends . 

The  procedure  employed  to  compute  the  variable  fre- 
quency magnitudes  was  described  in  detail  in  the  previous 
section  of  this  report.  As  described  there,  estimates  of 
magnitude  at  different  frequencies  are  based  on  the  peak 
amplitudes  of  the  envelopes  of  narrow  band  filters  applied 
to  the  signals  of  interest.  The  peak  signal  filter  ampli- 
tudes were  corrected  for  noise  based  on  estimates  made  from 
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Fiaure  3.1a.  Time  series  corresponding  to  an  mb  = 5.3  NTS  ex 
plosion.  Wagtail,  recorded  at  the  Yellowknife  array  in  Canada 


Time  (sec) 

Figure  3.1b.  YKA  short-period  recording  of  the  mb  = 
explosion,  Diana  Mist. 
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Figure  3.1c.  YKA  short-period  recording  of  an  i%  = 5. 
of  California  earthquake. 


3 Gulf 


Figure  3. Id.  YKA  short-period  recording  of  a small  (m^  = 4. 
Gulf  of  California  earthquake. 
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the  24  second  noise  windows  preceding  each  signal.  In  the 
present  VFM  discrimination  test  of  the  Yellowknife  data  base, 
twenty  filters  with  center  frequencies  ranging  from  0.25  Hz 
to  4.0  Hz  were  applied  to  each  noise  and  signal  time  series. 
The  preliminary  results  of  this  test  are  discussed  in  the 
following  sub-section  of  this  report. 


3 . 3 RESULTS  OF  THE  VFM  DISCRIMINATION  TEST 

As  mentioned  above,  each  of  the  signal  time  series 
was  filtered  by  twenty  narrow  band  filters  with  center  fre- 
quencies ranging  from  0.25  Hz  to  4.0  Hz.  The  frequency 
range  was  decided  on  after  examining  spectra  of  ten  different 
earthquakes  and  explosions.  Magnitude  estimates,  ir^(f), 
based  on  a narrow  band  filter  with  a relatively  low  center 

frequency,  f , were  then  plotted  versus  estimates  based  on 
c 

high  frequency  filters.  This  was  repeated  for  various  com- 
binations of  center  frequencies  and  all  of  the  plots  were 
then  examined  for  discrimination. 

Figure  3.2  shows  preliminary  results  for  the  North 
American  data  base.  The  filter  center  frequencies  that 
have  yielded  the  best  discrimination  results  to  date  are, 
as  indicated,  0.425  Hz  and  2.5  Hz.  In  this  figure  the  Gulf 
of  California  earthquakes  are  denoted  by  the  open  circles, 
the  NTS  explosions  by  the  x's.  In  general  most  of  the  earth- 
quakes on  this  figure,  which  are  reported  in  the  Preliminary 
Determination  of  Epiaenters  to  be  shallow  focus,  separate 
from  the  explosion  population.  However,  compared  with  the 
VFM  discrimination  results  previously  obtained  for  Eurasian 
events  recorded  at  LASA  (Fig.  3.3)  the  results  in  Fig.  3.2 
are  not  nearly  as  impressive  in  terms  of  complete  separation 
of  earthquake  and  explosion  populations. 

Before  we  can  draw  any  firm  conclusions  about  the 
data  set  in  Fig.  3.2  more  will  have  to  be  learned  about 
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Figure  3.3.  Spectral  magnitudes,  m^,  computed  at  0.45  Hz 

and  2.25  Hz.  The  presumed  explosions  numbered 
35  and  138  occurred  at  Novaya  Zemlya. 
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those  earthquakes  that  overlap  part  of  the  explosion  popula- 
tion. Note,  however,  that  a very  important  and  positive 
point  about  the  populations  i.r.  vhis  figure  is  the  complete 
separation  of  the  very  small  magnitude  events. 

3.4  SUMMARY 

A large  population  of  North  American  events,  namely 
Gulf  of  California  earthquakes  and  NTS  explosions,  were 
tested  for  discrimination  using  the  VFM  technique.  For 
certain  combinations  of  low  versus  high  filter  center  fre- 
quencies, a plot  of  the  variable  frequency  magnitude  esti- 
mates resulted  in  general  discrimination  of  the  earthquake 
and  explosion  populations,  with  some  overlap  of  the  two 
populations.  While  the  degree  of  discrimination  obtained 
to  date  is  less  for  the  North  American  event  populations  than 
for  Eurasian  event  populations  previously  studied,  it  must 
be  pointed  out  that  complete  separation  of  those  Gulf  of 
California  earthquakes  and  NTS  explosions  with  m^  (NOAA) 
values  between  4.0  and  4.5  was  achieved. 
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IV.  THE  EFFECT  OF  TECTONIC  RELEASE  ON  THE  TELESEISMIC 
SHORT  PERIOD  SEISMOGRAM 


4.1  INTRODUCTION 


Observations  of  the  seismic  waves  generated  by 
underground  nuclear  explosions  have  led  many  investigators 
to  conclude  that  the  explosion  is  often  accompanied  by  some 
type  of  induced  tectonic  stress  release.  An  excellent  sum- 
mary of  the  background  leading  to  this  conclusion  may  be 
found  in  Aki  and  Tsai  [1972],  There  are,  in  essence,  two 
explanations  of  the  predominant  mechanism  by  which  the  stress 
is  released.  The  first,  advanced  by  Archambeau  and  coworkers 
(Press  and  Archambeau  [1962],  Archambeau  [1964,  1968,  1972], 
Smith,  et  al.  [1969],  Lambert,  et  al . [1972],  Archambeau  and 
Sammis  [1970]),  assumes  that  the  mechanism  is  stress  relaxa- 
tion around  the  predominantly  spherical  fracture  zone  created 
by  the  explosion.  The  second  model  is  that  of  the  release  of 
stress  along  a pre-existing  fault  plane.  Evidence  for  the 
latter  model  is  summarized  by  Aki  and  Tsai  [1972]  where  it 
is  referred  to  us  the  "trigger  model".  In  either  case,  the 
center  of  dilatation  explosion  source  is  modified  by  the 
superposition  of  a double-couple  source  which  is  either  super- 
imposed on  the  explosion  (Archambeau1 s model)  or  somewhat 
removed  in  space  and  time. 

The  association  of  tectonic  stress  release  with  explo- 
sions is  important  for  a number  of  reasons.  Failure  to  ac- 
count for  its  influence  can  lead  to  bias  in  explosion  yield 
estimates  from  body  and  surface  wave  measurements.  Under 
certain  circumstances  it  is  conceivable  that  tectonic  re- 
lease could  cause  a failure  to  contain  the  radioactive  gases 
emitted  by  the  explosion.  Further,  an  accurate  estimate  of 
the  source  generated  elastic  wavelet  is  of  great  value  for 
determining  the  modulation  of  this  pulse  by  its  travel 
through  the  earth  and,  therefore,  attempts  to  invert  seismic 
data  to  uncover  the  structure  of  the  earth. 
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Interest  in  tectonic  release  has  been  primarily 
directed  to  an  explanation  of  the  commonly  observed  SH  and 
Love  waves  excited  by  explosions.  As  far  as  surface  waves 
are  concerned,  the  two  suggested  mechanisms  for  tectonic 
stress  release  are  quite  similar  and  both  have  been  found 
to  be  consistent  with  observations. 

The  possible  effect  of  tectonic  release  on  the  body 
wave  signature  of  explosions  is  also  of  interest.  Of  parti- 
cular significance  is  the  amplitude  from  which  body  wave 
magnitude  (m  ) is  determined.  It  is  our  intention  here  to 
identify  the  situations  in  which  tectonic  release  can  have  an 
important  effect  on  the  short  period  recording. 

Archambeau's  theory  has  been  cast  in  a form  suitable 
for  computation  with  the  controlling  parameters  being  re- 
lated to  the  properties  of  the  rock  at  the  hypocenter  and 
the  magnitude  and  orientation  of  the  prestress  field.  The 
computational  results  are  given  in  the  form  of  the  double- 
couple or  quadrupole  terms  in  the  multipolar  expansion  of 
the  displacement  field.  These  can  be  added  to  the  explosion 
generated  monopole  to  obtain  a complete  source  representation. 
Theoretical  seismograms  can  then  be  computed. 

Varying  the  parameters  in  Archambeau's  model,  we  can 
discover  the  values  which  must  obtain  in  order  for  tectonic 
release  to  have  an  effect  on  the  body  wave  recording.  Ap- 
plying other  evidence,  one  can  then  identify  those  situations 
in  which  these  parameters  may  be  reasonable  or  likely. 

The  double-couple  for  the  triggering  of  dislocation 
along  a nearby  fault  plane  is  mainly  different  (for  tele- 
seismic  body  waves)  from  Archambeau's  in  that  this  source  is 
separated  from  the  explosion  hypocenter.  This  quadrupole 
source  spectrum  would  be  expected  to  depend  on  stress  drop 
and  fault  (weak  zone)  dimensions  in  an  analogous  way.  While 
the  computations  to  be  shown  are  for  Archambeau's  model,  an 
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extrapolation  to  the  alternative  mechanism  is  straightforward. 
The  analysis  presented  is  therefore  relevant  to  the  general 
question  of  the  effect  of  tectonic  release  on  the  short  period 
P wave  record. 

4 . 2 ARCHAMBEAU 1 S TECTONIC  RELEASE  MODEL 

Before  proceeding,  a brief  discussion  of  Archambeau's 
model  and  the  important  parameters  controlling  the  size  of 
the  tectonic  release  component  of  the  radiation  field  is  in 
order.  The  physical  mechanism  on  which  Archambeau's  theory 
is  based  is  summarized  in  Fig.  4.1  (from  Archambeau  [1972]), 
The  zone  of  radius  Rc  is  made  up  of  pulverized  material 
and  is  assumed  to  be  created  at  a rate  denoted  V . The 

I\  j 

radially  cracked  zone  of  radius  R^  is  assumed  to  be  created 

at  a rate  V . The  zone  from  which  the  tectonic  prestress 
R 

is  released  is  of  radius  Rg.  A pure  shear  prestress  field 
of  amplitude  o{0>  is  applied. 

Due  to  the  filtering  properties  of  the  earth  and  the 
short  period  instrumentation,  teleseismic  body  wave  observa- 
tions generally  include  only  a narrow  sampling  of  the  source 
spectrum,  between  about  0.5  and  2.0  Hz.  The  important  para- 
meters for  the  tectonic  release  are  then  ^ and  R^  as 

the  amplitude  of  the  spectrum  is  directly  proportional  to 

a<°>*  and  R3 . Note  that  since  the  shatter  zone  volume 
0 , 
should  be  proportional  to  yield  and  to  R#,  the  tectonic 

release  amplitude  is  directly  proportional  to  yield.  The 

corner  frequency  depends  on  V_  and  R but  is  generally 

2 0 

sufficiently  high  to  be  of  little  importance.  Archambeau's 
assumption  of  a finite  R , which  is  interpreted  as  represent- 
ing  a nonuniform  prestress  field,  has  led  to  considerable 

*The  proportionality  is  more  fundamentally  to  the  strain 
which  equals  the  quotient  of  the  stress  and  the  material 
rigidity. 
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Figure  4. 
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. Schematic  representation  of  the  radial  zone  of 
material  around  a large  explosion  in  a pre- 
stressed medium.  Rc  and  Rq  are  radii  enclosing 
zones  of  failure  and  vary  with  explosive  yield, 
prestress  and  medium  type  while  Rg  is  the 
radius  within  which  most  of  the  relaxation  ef- 
fects occur  (from  Archambeau  [1972] ) . 
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controversy  (sec,  for  example,  Snoke  [1975] ) since  the  re- 
sulting source  spectrum  is  peaked.  While  the  assumed  value 

of  R is  of  paramount  importance  for  SH  and  Love  wave 
s 

studies  it  turns  out  to  be  of  little  consequence  to  the 
study  undertaken  here. 

In  subsequent  sections  we  will  first  describe  the 
character  of  the  source  radiation  for  an  explosion  plus 
tectonic  release.  Synthetic  seismograms  will  then  be  com- 
puted and  studied  to  determine  the  effect  of  the  tectonic 
release. 

4.3  EXPLOSION  MODELING 

Our  objective  is  to  determine  the  effect  of  the  addi- 
tion of  a tectonic  release  generated  quadrupole  component  to 
the  monopole  explosion.  The  amplitude  of  the  elastic  waves 
generated  by  the  explosion  itself  is  then  obviously  of 
critical  importance.  This  portion  of  the  source  is  discus- 
sed in  this  section. 

The  equivalent  elastic  source  for  the  spherically  sym- 
metric explosion  is  conveniently  discussed  in  terms  of  the 
reduced  displacement  potential  which,  in  the  frequency  domain, 
is  related  to  the  monopole  by 


Ay'  (u)  = -k* 


2 V(m) 
P vD 


(4.1) 


It  should  be  pointed  out  that  the  far  field  component  of  the 
displacement  spectrum  (in  a homogeneous  whole  space)  is  re- 
lated to  the  reduced  displacement  potential  by 


-if** 

P 


— ikpR 


(4.2) 


Deterministic  computational  methods  for  determining 
the  reduced  displacement  potential  for  underground  nuclear 
explosions  have  been  developed  by  J.  T.  Cherry  and  his  col- 
leagues (Cherry,  et_  al.  [1973,  1974al,  Bache,  et  al.  [1975b]). 
These  computations  employ  finite  difference  methods  and  re- 
quire estimates  of  the  constitutive  behavior  of  the  rock 
over  tile  range  of  stresses  encountered.  The  material  be- 
havior accounted  for  includes  that  attributable  to  water 
content,  air-filled  porosity,  brittle  failure  and  plastic 
flow.  Computed  source  representations  for  various  rock  types 
encountered  at  NTS  have  been  shown  to  be  consistent  with  both 
near-source  and  teleseismic  body  wave  data  (Bache,  et  al. 

[1975b] ) . 

/s 

In  Fig.  4.2  is  shown  the  amplitude  of  H' (w)  for  two 
calculations,  both  representing  NTS  tuff.  The  second  was 
computed  specifically  for  a class  of  explosions  in  saturated 
tuff  at  Yucca  Flat,  NTS.  The  peaking  in  the  source  spectrum 
for  the  latter  is  primarily  due  to  the  form  of  the  effective 
stress  law  in  the  constitutive  model  (see  Bache,  et  al.  [1975b] 
for  details) . 

The  source  spectra  of  Fig.  4.2  shall  be  used  to  repre- 
sent the  explosion  in  the  calculations  to  follow.  The  range 
of  | ’F  (oj)  | expected  to  be  encountered  at  NTS  is  indicated  in 
Fig.  4.3  which  is  taken  from  Bache,  et  al.  [1975a],  Note  that 
in  the  teleseismic  frequency  band,  Source  A is  considerably 
smaller  in  amplitude  than  Source  B. 

Since  the  calculations  leading  to  the  source  spectra 
of  Fig.  4.3  were  consistent  with  the  available  near  field 
data  as  well  as  teleseismic  body  wave  observations,  they  should 
provide  useful  guidance  in  selecting  the  important  parameter, 
R0'  for  tlle  tectonic  release.  In  particular,  the  greatest 
radius  at  which  material  yielding  occurred  was  monitored.  This 
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Source  B of  Fig.  4.2. 
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number,  denoted  R0,  should  provide  an  upper  limit  to  R 
The  computed  values  of  R0  are  summarized  in  Table  4.1. 


TABLE  4.1.  COMPUTED  MAXIMUM  RADIUS  AT  WHICH  YIELDING  OCCURS 
FOR  150  KT 


Source 


Radius  (meters) 


A 

B (Yucca  Flat  Wet  Tuff) 
Yucca  Flat  Dry  Tuff 
Area  12 

Pahute  Mesa  Rhyolite 


490 

670 

310 

850 

360 


4.4  THE  EQUIVALENT  ELASTIC  SOURCE 

A compatible  tectonic  release  quadrupole  component  is 
now  added  to  the  explosion  monopole  and  the  resulting  elastic 
field  is  examined.  The  specific  example  chosen  is  a 150  kt 
explosion  in  tuff.  The  important  features  of  the  four  source 
representations  to  be  studied  are  summarized  in  Table  4.2. 

Tuff  is,  of  course,  a common  material  for  underground 
tests  at  NTS.  The  yield  of  150  kt  is  typical  of  many  such 
tests.  While  we  will  be  dealing  with  this  specific  case, 
the  results  will  be  representative  of  a wide  range  of  tests 
and  near  source  materials,  as  will  be  pointed  out  in  the  en- 
suing discussion. 

The  selection  of  parameters  for  the  tectonic  release 
calculations  was  based  on  the  following  considerations.  The 
measured  cavity  radius  for  150  kt  shots  in  tuff  is  on  the 
order  of  60  meters.  Recalling  that  the  model  is  rather  in- 
sensitive to  Rc,  a value  of  60  meters  was  then  used  in  all 
calculations.  As  mentioned  in  the  previous  section,  an 
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estimate  of  R is  given  by  the  computed  outer  radius  at 
o 

which  yielding  occurs.  Since  the  tectonic  release  scales 

with  R , we  shall  be  interested  in  a range  of  values  and 

o 

particularly  the  minimum  R#  required  for  tectonic  release 
to  affect  the  body  wave.  The  same  is  true  for  the  prestress 
which  is  initially  chosen  to  correspond  to  a strain  of 
» 10“ \ a value  not  unreasonable  for  NTS  materials.  Note 
that  the  tectonic  release  component  for  Sources  1 and  4 is 
essentially  identical  since  it  is  prestrain  which  is  the 
determining  parameter. 

The  elastic  radiation  field  generated  by  the  four 
sources  of  Table  4.2  will  now  be  studied  by  examining  far  field 
spectra  and  radiation  patterns.  The  far  field  component  of 
the  radiation  pattern  is  that  which  decays  as  R . The  co- 
ordinate system  to  be  used  is  depicted  in  Fig.  4.4. 

First,  consider  the  case  of  a pure  strike-slip  tec- 
tonic release.  That  is,  let  o =*  a(0).  Spectra  for  Source* 

xy 

1,  2,  and  3 at  azimuth  <t>  = 45°  and  two  takeoff  angles, 
t « 15°,  45°,  are  shown  in  Fig.  4.5.  The  radiation  field  is 
composed  of  P,  SV  and  SH  components  (the  outgoing  S wave  is 
polarized  with  respect  to  the  free  surface) . For  Source  1 
the  1.0  Hz  spectral  component  radiation  patterns  are  plotted 
in  Fig.  4.6  at  azimuth  45°  and  takeoff  angle  15°.  The  azimuth 
45°  is  a maximum  on  the  SV  radiation  pattern  and  a node  for 
SH  waves.  Since  we  will  not  be  concerned  with  SH  and  Love 
waves,  the  SH  component  will  not  be  shown.  Finally,  in 
Fig.  4.7  the  1.0  Hz  radiation  patterns  for  Source  2 are  shown. 

For  teleseismic  body  waves  the  important  energy  leaving 
the  source  is  at  takeoff  angles  = 5°-20°  for  the  direct  P wave 
and  * 160°-175°  for  the  free  surface  reflected  pP  and  sP 
phases.  For  a takeoff  angle  of  15°,  the  ray  parameter 
p a 0.11  sec/km  and  the  free  surface  reflection  coefficients 
are  approximately  -0.92  and  0.16  for  pP  and  sP,  respectively 
(Cerveny  and  Ravindra  (1971),  Section  2.4). 
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Figure  4.7.  Normalized  radiation  patterns  for  the  1.0  Hz  spectral  component  for 
Source  2.  The  normalization  is  as  in  the  previous  figure. 
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Concerning  Figs.  4. 5-4. 7,  the  following  observations 
are  made: 

• For  the  spectra  shown,  the  SV  wave  is  about  the 

same  size  as  the  P wave  in  the  band  0.5  < f < 2.0  Hz. 
It  must  be  much  larger  for  sP  to  strongly  affect  the 
body  wave  record. 

• At  the  relevant  takeoff  angles,  the  P wave  is 
little  affected  by  tectonic  release. 

• The  portion  of  the  spectrum  which  is  affected  by 

Rg  is  outside  the  band  of  interest  for  body  waves. 

However,  R is  clearly  important  for  long  period 
s 

waves. 

• The  effect  of  R is  primarily  to  scale  the 
spectrum  with  r|.  The  spectrum  is  also  shifted 
to  lower  frequencies  with  increasing  R#,  but 
this  effect  is  minor  for  the  range  of  appropriate 
values. 

• The  strike-slip  orientation  is  rather  unfavorable 
for  observing  tectonic  release  in  teleseismic 
body  waves.  Other  orientations  of  the  prestress 
field  are  considerably  more  favorable. 

The  tectonic  release  controlling  parameters  R 

(0)  0 
and  a used  above  are  large  but  seem  not  unreasonable 

for  NTS  tuff  explosions.  We  should  also  attempt  to  define 

a source  which  includes  the  greatest  values  that  can  be 

justified  for  this  model.  Therefore,  let  us  multiply  the 

tectonic  release  component  in  Source  1 by  four.  Since  the 

(O' 

amplitude  is  essentially  directly  proportional  to  a ’ and 

R3,  this  is  equivalent  to  taking  = 240  bars,  R = 0.5  km 

0 (0)  0 
or  a =60  bars,  R^  = 0.79  km  or  various  other  combinations 

giving  a ^ R3  30.  The  Source  1 spectrum  at  $ = 45°, 
o 

t = 15°,  is  compared  to  that  for  the  same  source  with  the 
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tectonic  release  multiplied  by  four  in  Fig.  4.8.  Parenthe- 
tically it  should  be  noted  that  it  is  difficult,  but  not 

impossible,  to  justify  using  values  of  R*  > 30  in  the 

o 

calculations.  On  the  other  hand,  one  might  argue  that 
0 ^ R30  must  be  less  than  1.  The  physics  is  rather  poorly 

understood  and  remains  a subject  for  conjecture.  Studies 
of  the  effect  of  tectonic  release  provide  the  best  hope  for 
fixing  these  parameters. 

4.5  COMPUTATION  OF  THEORETICAL  SEISMOGRAMS 

Having  developed  some  understanding  of  the  source 
displacement  spectra,  we  now  shift  attention  to  the  compu- 
tation of  theoretical  seismograms  to  see  directly  the  effect 
of  tectonic  release.  The  computations  are  based  on  the 
representation  of  an  equivalent  elastic  source  buried  in  a 
stack  of  plane  elastic  layers,  which  was  developed  in  Ap- 
pendix B,  Bache,  et:  al . [1975a] . Input  to  the  representation 
includes  the  crustal  structure,  depth  of  burial,  the  multi- 
polar coefficients  defining  the  source  and  the  ray  parameter 
of  interest  for  the  teleseismic  wave.  The  output  is  the 
displacement  spectrum  along  the  specified  ray. 

We  are  interested  in  isolating  the  source,  which  is 
taken  to  be  the  explosion  plus  tectonic  release,  and  its 
modulation  by  the  near  source  geology.  Therefore,  while  we 
need  realistic  seismograms,  modification  by  the  upper  mantle 
and  receiver  region  crustal  reverberations  are  minimized  in 
the  computations.  The  following  elements  are  included: 

• A constant,  G , is  used  to  represent  the  geo- 

. 5 

metric  spreading  in  the  upper  mantle.  This 
constant  is  computed  using  the  method  of  Julian 
and  Anderson  [1968]  and  is  dependent  on  the 
upper  mantle  velocity  model  and  epicentral 

o 

distance  chosen.  The  earth  model  CIT109 
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Figure  4. 


. Normalized  displacement  spectra  at  <J>  = 45°, 
x = 15°  for  Source  1 compared  to  the  spectra 
for  the  same  source  with  the  tectonic  release 
component  quadrupled. 


Here  f is  frequency  and  the  controlling 
parameter  T/Q  is  the  ratio  of  travel  time 
(T)  to  the  average  path  material  quality 
factor  (Q) . For  all  computations  T/Q  = 0.9. 


t Crustal  reverberations  at  the  receiver  are 
computed  using  the  Thompson-Haskell  method 
(e.g.,  Haskell  [1962]).  A crustal  model 
having  no  very  strong  amplitude  effects  was 
chosen  fo’'  all  computations.  The  model  is 
shown  in  Fig.  4.9. 

• The  transfer  function  for  a standard  short 
period  instrument,  the  LRSM  Benioff,  was 
applied. 

As  the  illustrative  example,  we  have  chosen  a 150  kt 
explosion  in  NTS  tuff.  For  all  calculations  the  depth  of 
burial  was  taken  to  be  0.8  km,  a value  typical  of  explosions 
in  this  yield  range.  Two  crustal  models  were  used  for  the 
source  region  and  these  are  shown  in  Fig.  4.9.  The  first 
model  (SSI)  has  a very  strong  impedance  contrast  between  the 
tuff  and  underlying  hard  rock.  A tuf f-paleozoic  rock  con- 
trast of  this  kind  is  found  at  Yucca  Flat.  The  second  model 
(SS2)  has  this  contrast  spread  over  a wavelength  or  so  (at 
2 1 HZ)  . 
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Figure  4. 


P Wave  Velocity  (km/sec) 


1.0  2.0  3.0  4.0  5.0  6.0 


. Compressional  wave  velocity  versus  depth  for 
three  crustal  models  used  in  the  computations 
the  source  region  models  SSI  and  SS2  and  the 
receiver  region  model  RS. 
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Synthetic  seismograms  for  the  explosion  alone  (no 
tectonic  release)  are  shown  in  Fig.  4.10.  As  an  indication 
that  these  theoretical  records  are  representative  of  reality, 
several  Palmer  Network  observations  of  the  Starwort  event 
are  also  shown  in  the  figure.  Starwort  was  an  85  kt  shot, 
buried  at  0.564  km  in  tuff  at  Yucca  Flat.  Taking  account 
of  the  shallower  burial  depth,  we  see  that  the  frequency 
content  and  shape  of  the  first  few  swings  on  the  records  are 

comparable. 

In  Figs.  4.11  and  4.12  are  shown  the  effect  of  in- 
cluding the  strike-slip  tectonic  release  component  discussed 
in  the  previous  section.  The  tectonic  release  model  is  that 
designated  No.  1 in  Table  4.2  and  the  same  model  with  the 
tectonic  release  component  multiplied  by  two  and  by  four. 

The  latter  is  taken  to  represent  the  maximum  credible  amount 
of  tectonic  release  for  the  physical  model  considered.  Com- 
putations were  made  with  both  the  crustal  models  SSI  and 
SS2.  In  view  of  the  radiation  patterns  (Fig.  4.5),  only  the 
azimuths  45°  and  135°  were  considered  since  the  tectonic 
release  influence  is  maximized  at  tiiese  azimuths. 

It  is  immediately  clear  from  Figs.  4.11  and  4.12 
that  a pun  strike-slip  toctonic  release  has  insignificant 
influence  on  teleseismic  body  waves.  Even  for  the  maximum 
amount  (A  = 4)  at  the  optimal  azimuths,  the  effect  is  hardly 
noticeable,  especially  in  view  of  the  variation  between  ob- 
servations at  nearby  stations  seen  in  Fig.  4.10. 

4.6  TECTONIC  RELEASE  FROM  AN  OPTIMALLY  ORIENTED  PRESTRESS 

FIELD 

It  is  natural  to  inquire,  what  is  the  effect  of  tec- 
tonic release  that  is  oriented  to  maximize  the  effect  on 
teleseismic  short  period  observations?  Such  an  orientation 
results  from  a pure  shear  prestress  field  that  dips  at  45°. 


‘ 


(b)  Observed  seismograms  fran  the  Palmar  Network  in  Alaska  for  the 
Yucca  Flat  event  Starwort.  The  station  codes  are  indicated  as 
is  the  magnification  at  1 Hz. 

Figure  4.10.  Synthetic  seismogr^ns  fo^  the  explosion  alone  (no 
tectonic  release)  and  Palmer  Network  observations 
of  an  event  similar  to  that  modeled. 
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In  the  coordinate  system  of  Fig.  4.4,  one  such  field  is 
given  by  a^z  = "°Xy  “ //I.  This  corresponds  to  rotating 

the  strike-slip  sources  of  the  previous  section  by  45°  about 
the  y-axis.  Sample  radiation  patterns  for  the  new  source  are 
shown  in  Fig.  4.13.  For  this  source  the  optimal  azimuths  for 
tectonic  release  effects  should  be  + 30°  from  the  x-axis. 

Theoretical  seismograms  including  tectonic  release 
at  the  "optimal"  orientation  are  shown  in  Figs.  4.14  and 
4.15.  Once  again  we  have  used  Source  1 of  Table  4.2  and  the 
same  source  *.ith  the  tectonic  release  component  doubled  and 
quadrupled.  The  following  conclusions  may  be  drawn: 

9 For  the  P.  =0.5  km,  cr^  = 60  bars  source, 
o 

the  influence  of  tectonic  release  is  rather 
insignificant,  even  at  these  optimal  azimuths 
for  the  optimal  orientation. 

• For  a doubled  or  quadrupled  tectonic  release, 
the  effects  become  important.  In  particular, 
a strong  variation  with  azimuth  becomes  ap- 
parent. The  maximum  amplitude  from  which 
would  be  measured  varies  by  as  much  as  a factor 
of  two,  dependent  on  azimuth. 

So  far  we  have  restricted  attention  to  Sources  1-3  of 
Table  4.2  and  thus  the  explosion  model  denoted  A in  Fig.  4.2., 
To  be  certain  that  our  conclusions  have  not  been  unduly  in- 
fluenced by  this  choice,  we  now  consider  Source  4 of  Table  4.2, 
which  is  quite  similar  to  Source  1 except  that  the  explosion 
model  B is  used.  The  relevant  synthetic  seismograms  are 
shown  in  Fig.  4.16.  For  this  example  we  only  show  the  cases 
where  the  effect  of  tectonic  release  is  maximized.  We  find 
that  it  is  only  at  the  optimum  azimuths  for  the  optimum 
orientation,  that  tectonic  release  has  a noticeable  effect. 

This  is  expected  since  the  explosion  induced  elastic  wave 
(see  Fig.  4.2)  is  relatively  larger  than  that  for  Source  1 
while  the  tectonic  release  component  remains  the  same. 
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This  theoretical  analysis  represents  an  attempt  to 
answer  the  question;  under  what  circumstances  can  tectonic 
release  have  an  important  effect  on  the  telesei9mic  body 
wave  signature  of  explosions?  Underlying  the  analysis  are 
several  fundamental  assumptions:  (1)  The  source  is  composed 

of  a spherically  symmetric  radiator  of  elastic  waves,  the 
explosion,  plus  a double-couple  representing  tectonic  re- 
lease. (2)  The  computed  reduced  displacement  potential 
representations  of  the  explosion  (Figs.  4.2  and  4.3)  are 
reasonably  accurate.  (3)  The  earth  models  and  computational 
techniques  used  for  constructing  the  theoretical  seismograms 
are  sufficiently  realistic  to  allow  conclusions  tc  be  drawn 
about  the  contribution  of  the  tectonic  release  component  to 
the  teleseismic  record. 

With  the  assumptions  listed  above,  we  have  done  the 
following: 

• Since  the  important  factor  is  the  ratio  of  the 
explosion  component  of  the  displacement  field 
to  that  due  to  tectonic  release,  the  sample 
calculations  were  done  for  explosions  which 
couple  relatively  weakly  into  elastic  wmves 
(Figs.  4.2  and  4.3). 

• Using  Archambeau's  model  of  the  tectonic  stress 
release  into  the  explosion  created  shatter  zone, 
we  chose  what  seemed  to  be  the  largest  values 
that  could  be  justified  for  the  weak  zone  dimen- 
sions and  prestress  or  prestrain. 

• Two  orientations  of  the  prestress  were  invesgi- 
gated,  one  a pure  strike-slip  and  the  second 
rotated  45°  about  a horizontal  axis.  The  latter 
is  essentially  optimal  for  propagating  tectonic 
release  generated  body  waves  to  teleseismic 
distances . 
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• For  either  orientation,  theoretical  seismo- 
grams were  computed  at  azimuths  where  the 
tectonic  release  influence  was  strongest. 

Let  us  first  assume  that  the  mechanism  for  the  relaxa- 
tion of  tectonic  release  if  ti at  hypothesized  by  Archambeau. 
Then  we  can  conclude  that  the  body  wave  signature  of  explo- 
sions can  only  be  affected  if  the  prestress  field  is  oriented 
near  the  optimum.  That  is,  if  the  vertical  component  of  the 
prestress  is  of  the  same  order  as  the  horizontal  (e.g., 

| a j « |0  |).  Even  for  this  orientation  it  is  necessary 

that  the  tectonic  release  parameters  be  rather  large.  For 
the  150  kt  example  studied,  we  found  that  the  product  of 


the  prestress,  a'u/,  and  R3 , where  R is  the  radius  of 

6 0 

the  shatter  zone  within  which  the  material  rigidity 
essentially  vanishes  (y  = 0) , must  be  * 30  km-bar.  For 
explosions  characterized  by  a larger  reduced  displacement 
potential  (e.g.,  those  at  Pahute  Mesa),  these  parameters 
must  be  commensurately  larger.  For  such  cases  tectonic 
release  can  significantly  influence  the  shape  of  the  dis- 
placement pulse  character. zing  the  source  and  change  the 
amplitude  from  which  is  measured  by  a factor  of  two  or 
so.  These  effects  are,  of  course,  strongly  dependent  on 
azimuth.  Once  again,  it  is  only  at  optimal  azimuths  for 
the  optimal  orientation  and  the  most  favorable  ratio  of 
tectonic  release  to  explosion  generated  elastic  waves  that 
the  effect  becomes  important.  Therefore,  it  seems  reason- 
able to  dismiss  the  tectonic  release  mechanism  suggested 
by  Archambeau  as  an  important  contributor  to  the  body  wave 
signature  of  explosions. 

What  if  the  tectonic  release  mechanism  is  the  trigger 
model  proposed  by  Aki  and  Tsai?  First,  we  would  expect  the 
time  lag  between  the  explosion  and  tectonic  release  generated 
waves  to  be  greater  than  studied  here.  Otherwise,  the  source 
characteristics  important  for  tcieseismic  body  waves  should 
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be  reasonably  modeled  by  Archambeau's  source  theory.  In  fact, 
the  double-couple  term  for  the  explosion  induced  tectonic  re- 
lease is  essentially  identical  to  the  leading  term  in  the 
multipolar  expansion  for  the  earthquake  source  theory  de- 
veloped by  Archambeau  and  Minster  (Archambeau  [1964] , 

Minster  [1973] ) . Here  is  taken  to  be  representative 

of  half  the  fault  length  and  a(0)  is  the  stress  drop. 

From  the  theoretical  seismograms  for  the  weak  coupling 
150  kt  shot  studied,  we  know  that  tectonic  release  has  almost 
no  effect  on  the  short  period  P wave  record  when  o^^R3  » 7,5. 
no  mattor  what  the  orientation.  For  optimal  orientations 
significant  effects  could  be  noticed  when  o^^R3  z 30.  Aki 
and  Tsai  [1972]  conclude  that  the  tectonic  release  associated 
with  large  NTS  explosions  is  a low  stress  drop  (s  10  bars) 
phenomenon.  For  o(0)R3  = 30,  we  then  would  expect  a signi- 
ficant contribution  to  the  body  waves  for  fault  lengths  on 
the  order  of  3 km.  Aki  and  Tsai  postulate  fault  lengths  on 
the  order  of  4 km  for  the  235  kt  shot  Bilby  and  10  km  for 
Boxcar  (1200  kt)  and  Benham  (1100  ktr  data  from  Springer  and 
Kinnaman  [1971]).  Scaled  to  150  kt  and  halved,  these  are 
equivalent  to  R#  » 1.7  km  for  Bilby  and  s 2.6  km  for  the  two 
larger  events.  Then  if  the  stress  drop  is  20-30  bars  ov<ir 
faults  of  this  order  length,  and  if  the  fault  displacement 
is  neither  strike-slip  nor  dip-slip  but  somewhere  between  and 
is  associated  wi  h low  coupling  explosions,  we  would  expect 
a significant  contribution  to  the  short  periol  P wave  record. 
Release  of  tect  '.lie  stress  by  this  mechanism  is  therefore 
somewhat  more  likely  to  contribute  to  the  short  period  record 
than  by  the  mechanism  proposed  by  Archambeau.  Still,  a 
coincidence  of  a number  of  favorable  circumstances  is  re- 
quired and  this  seems  unlikely  for  moot  events. 

With  this  analysis  we  are  forced  to  conclude  that  for 
most,  perhaps  almost  all,  events  .it  is  safe  to  ignore  tec- 
tonic release  as  a contributor  to  the  short  period  P wave 
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recording.  Still,  a number  of  questions  remain  that  can  only 
be  answered  with  further  study.  The  choice  between  the  trig- 
ger model  and  Archambeau's  model  of  tectonic  release  may  be 
resolved  by  a study  of  short  period  SH  waves  since  depth  of 
burial  should  be  diagnostic.  The  amplitude  of  the  tectonic 
release  component  could  also  be  inferred  from  such  a study. 

The  amplitude  of  the  tectonic  release  component  is 
best  determined  from  surface  wave  observations.  The  work 
of  Toksoz  and  Kehrer  [1972]  represents  a good  beginning. 

Recent  developments  in  the  computation  of  the  explosion 
generated  portion  of  the  source  (Bache,  et  al.  [1975b]) 
should  be  quite  helpful  in  further  study. 

Finally,  it  shou  . be  mentioned  that  source  asymmetries, 
whether  due  to  emplacement  configuration  or  working  point 
media  inhomogeneities  (e.g.,  jointing),  can  also  cause  the 
superposition  of  a double-couple  on  the  source  monopole  (for 
example,  see  Cherry,  et  aL  [1974]).  The  effect  is  virtually 
indistinguishable  from  that  due  to  Archambeau's  model  of  tec- 
tonic release.  It  may  be  that  all  three  mechanisms  for 
generating  the  double-couple  are  present  to  some  degree  with 
most  explosions. 
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V.  A THREE-DIMENSIONAL  FINITE  DIFFERENCE  SIMULATION  OF 
STICK-SLIP  EARTHQUAKE  FAULTING  - THE  EQUIVALENT 

ELASTIC  SOURCE 

5.1  INTRODUCTION 

One  of  the  important  objectives  of  this  contract  is 
to  uncover  the  nature  of  the  earthquake  source,  especially 
those  features  that  discriminate  between  earthquakes  and 
underground  explosions.  With  this  objective  in  mind,  a 
stick-slip  rupture  model  has  been  incorporated  into  a three- 
dimensional  Lagrangian  finite  difference  code.  The  model 
has  been  exercised  to  compute  the  radiated  energy  due  to  a 
particular  fault  geometry.  In  this  case  the  geometry  was 
that  of  a bilateral  rupture  along  a 1 km  fault  length. 

A powerful  technique  for  studying  the  elastic  waves 
radiated  from  nonlinear  source  regions  is  first  to  replace 
the  source  with  an  equivalent  elastic  point  source.  The 
computational  techniques  of  theoretical  seismology  can  then 


be  applied.  The  computation  of  the  equivalent  elastic  source 
from  finite  difference  source  calculations  was  described  in 
Appendix  G of  our  last  report  under  this  contract  (Bache, 
et  al^  [1975a] ) . The  point  source  representation  is  unique 
and  exact  in  the  sense  that  everywhere  outside  the  "elastic 
radius"  at  which  the  point  source  was  obtained,  the  radiation 
from  the  equivalent  point  source  is  identical  to  that  from 
the  original  nonlinear  source. 

In  this  section  the  three-dimensional  finite  difference 
earthquake  (FDEQ)  calculation  is  described  with  emphasis  on  the 
equivalent  elastic  source  for  this  earthquake.  Of  particular 
interest  is  the  description  in  Section  5.5  of  the  far  field 
(propagating  to  teleseismic  distances)  displacement  field 
radiated  by  the  earthquake. 
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The  important  issue  is  what  this  finite  difference 
earthquake  simulation  tells  us  about  the  nature  of  the 
earthquake  source.  To  this  end  we  are  comparing  the 
character  of  this  source  to  data  and  to  analytical  models 
developed  by  others.  A number  of  interesting  results  are 
being  obtained  and  these  will  be  discussed  in  forthcoming 
reports. 


5.2  FAULT  CONFIGURATION  FOR  THE  THREE-DIMENSIONAL  CAL- 
CULATION 

A stick-slip  rupture  model  has  been  incorporated  into 
both  two-dimensional  (plane  strain)  and  three-dimensional 
Lagrangian  stress  wave  codes  by  Cherry  and  co-workers.  Re- 
sults of  calculations  of  near-field  earthquake  ground  motion 
from  the  plane  strain  version  of  the  model  were  reported  in 
Cherry  [1973] . As  reported  by  Cherry,  et  al.  [1374b] , this 
version  of  the  model  has  been  exercised  to  uncover  the  depen- 
dence of  peak  ground  motion  and  response  spectra  on  fault 
length,  rupture  velocity,  and  dynamic  stress  drop  during  rup- 
ture. Appendix  A of  Cherry  [1974]  shows  how  the  fault  sur- 
face is  isolated  in  terms  of  the  normal  and  tangential  compo- 
nents of  stress  occurring  at  the  interface.  In  Appendix  B 
of  the  same  reference,  contact  discontinuity  boundary  condi- 
tions are  applied  at  the  interface  and  the  difference  equa- 
tions required  to  simulate  the  fault  surface  during  rupture 
are  given. 

The  basic  mechanism  for  releasing  the  strain  energy 
in  this  rupture  model  is  the  relaxation  of  the  tangential 
stress  at  the  slipping  interface  (fault  surface)  from  its 
welded  value  to  its  kinetic  friction  value.  This  relaxation 
allows  adjacent  points  on  the  interface  to  nove  apart  (slip) 
The  rupture  heals  (adjacent  points  on  the  interface  stick) , 
if  the  relative  velocity  between  two  adjacent  points  change 
sign  and  if  the  tangential  stress  at  the  interface  is  suf- 
ficient to  maintain  continuity  of  tangential  velocity. 
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In  the  2-D  earthquake  model  the  fault  plane  is  of 
infinite  extent  in  the  direction  perpendicular  to  the  rupture 
direction  while  in  the  3-D  calculations  the  fault  plane  is  of 
finite  (specified)  extent.  In  comparing  calculations  done 
to  date,  there  is  another  significant  difference.  For  all 
the  2-D  calculations  the  faulting  was  unilateral  while  for 
the  single  3-D  calculation  done  to  date  the  rupture  was  bi- 
lateral. This  three-dimensional  simulation  of  bilateral 
earthquake  faulting  is  of  primary  interest  here. 

The  configuration  of  the  3-D  rupture  model  is  illus- 
trated schematically  in  Fig.  5.1.  The  primary  model  para- 
meters characterizing  the  3-D  stick-slip  earthquake  source 
calculation  are  summarized  below. 

• The  fault  is  located  at  the  origin  of  a 
Cartesian  (rectangular)  coordinate  system 
as  shown  in  Fig.  5.1. 

• The  fault  plane  for  the  rupture  model  is  in 
the  X-Y  plane. 

• The  rupture  velocity  has  an  average  magnitude 
of  VR  = 2.7  km/sec  and  is  parallel  to  the  Y- 
axis  and  in  the  positive  Y direction  for  Y > 0 
and  in  the  negative  Y direction  for  Y < 0,  cor- 
responding to  a bilateral  rupture.  The  direction 
of  particle  motion  is  indicated  in  the  figure. 

• The  fault  zone  is  defined  by  -0.5  <_  Y < 0.5 
and  -0.3  £ X £ 0.3,  in  units  of  kilometers, 
corresponding  to  a fault  of  length  (L)  1.0  km 
and  width  (W)  0.6  km. 

• The  prestress  is  composed  of  a single  com- 
ponent, ®2Y  = However,  the  stress 

drop  across  the  fault  boundary  is  taken  to 
be  0.15  kbar. 


Figure  5.1.  Configuration  of  the  bilateral  rupture  model 
with  fault  length  L,  width  Wf  and  prestress 
component  Sgy  as  used  in  the  3-D  finite  dif- 
ference earthquake  source  calculation.  The 
relation  between  the  Cartesian  (X,Y,Z)  and 
spherical  (R,0,<J>)  source  coordinate  systems 
is  illustrated. 
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The  3-D  finite  difference  calculation  used  a 
rectangular  grid  system  extending  from  0.0  to 
6.472  km  in  the  X and  Y directions  and  from 
-6.427  to  +6.427  km  in  the  Z direction.  A 
total  of  50  x 50  * 100  = 250,000  cells  were 
required.  The  rectangular  grid  was  evenly 
spaced  with  a grid  length  of  0.1  km  -.rom 
0.0  to  3.0  km  in  all  three  directions;  the 
grid  was  progressively  stretched  beyone  3.0  km. 

The  time  step  was  fixed  at  0.005  sec  and  the 
numerical  source  calculation  was  run  out  to  a 
total  elapsed  time  of  1.2150  sec  or  a total 
of  243  cycles. 

A total  of  319  monitoring  stations  were  arranged 
in  the  neighborhood  of  the  spherical  surface  of 
radius  R = 1.5  km  extending  over  one  octant  de- 
fined by  0°  < 6 ± 90°  and  0°  <_  <|>  £ 90°,  where 
6 is  the  colatitude  measured  from  the  Z-axis 
downward  to  a point  on  the  sphere  and  4>  is  the 
azimuthal  angle  measured  from  the  XZ  plane. 

The  primary  source  variables,  namely  the  dila- 
tation and  rotation  displacement  potentials, 
at  the  centers  of  the  eight  rectangular  cells 
surrounding  each  station  were  sampled  for  every 
time  step.  Also,  the  particle  displacement  and 
velocity  field  components  were  sampled  at  each 
monitoring  station. 

> The  material  properties  in  the  elastic  region  are 


P-wave  velocity;  vp  = 5.7  km/sec 
S-wave  velocity*  vg  = 3.4  km/sec 


Density: 

Bulk  modulus; 
Shear  modulus; 


p = 2.8  g/cc 
K = 478  kbar 
]i  = 324  kbar 
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As  a consequence  of  the  symmetry  properties  of  the 
3-D  stick-slip  earthquake  source,  a set  of  symmetry  relations 
can  be  derived  for  the  dilatation  and  rotation  potentials. 

By  applying  these  symmetry  properties,  numerical  values  of 
the  potentials  in  the  2nd  - 8th  octants  on  the  spherical 
surface  can  be  related  to  values  on  the  first  octant. 
Therefore,  the  set  of  monitoring  stations  can  be  restricted 
to  a single  octant  in  three-dimensional  space.  This  re- 
striction allows  a substantial  reduction  in  the  total  com- 
puting expenditure  required  to  carry  out  a full  3-D 
deterministic  simulation  of  a stick-slip  earthquake  without 
any  compromise  of  physical  reality. 

The  symmetry  properties  and  related  displacement 
boundary  conditions  for  the  stick-slip  bilateral  rupture 
model  are  discussed  in  Appendix  A.  The  derivation  of  the 
symmetry  properties  of  the  displacement  potentials  is  also 
given  in  the  same  Appendix,  along  with  a discussion  of  the 
special  computational  considerations  required  to  treat 
stations  near  the  boundary  surfaces,  represented  by  the  XZ, 
YZ,  and  XY  planes.  Finally,  the  symmetry  properties  re- 
lating values  of  both  displacement  and  potential  fields 
to  those  in  the  first  octant  are  summarized  in  Tables  A.l  ~ 
A. 3. 


5.3  EQUIVALENT  ELASTIC  SOURCE  REPRESENTATION 

Fundamental  to  studies  of  the  far-field  signature  of 
the  3-D  finite  difference  simulation  of  a stick-slip  earth- 
quake is  a representation  of  the  rupture  model  in  terms  of 
an  equivalent  elastic  source.  The  representation  of  a 
spatially  limited  source  in  terns  of  a multipole  expansion 
in  spherical  harmonics  has  been  discussed  in  a number  of  S* 
reports  (e.g.,  Cherry,  et  al.  [1975];  Bache,  et  ah  [1974] 
and  Bache,  et  al . [1975a]).  Only  the  final  analytic 
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expressions  will  be  introduced  here  to  provide  a convenient 
starting  point  and  to  establish  the  notational  conventions 
for  the  following  discussions. 

The  Fourier  transform  of  the  displacement  potentials 
is  expanded  in  terms  of  spherical  eigensolutions  of  the  wave 
equation  in  a homogeneous,  isotropic  linearly  elastic  medium 
as  follows: 


« 1 1 


(2) 


la  (r,0,<f>,u)i  - ^ ^ h^  ' (kar) 


£=0  m=0  s=0 


<“>  (e'*>  • 


(5.1) 


where  are  multipole  coefficient,  h|2^  are  spherical 

Hankel  functions  of  the  second  kind,  and  Y^ms  are  (unnormalized) 
spherical  harmonics  (Mo-rse  and  Feshbach  [1953]),  namely 


P™  (cosQ)  co.  mi}),  s=0] 


YMs  (e'^  = 


P?  (COS0)  sin  mi}),  s=ll 


(5.2) 


where  P^  are  associated  Legendre  functions.  The  multipole 
coefficients  in  the  time  domain  can  be  expressed  in  terms  of 
the  potentials  as 
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where  the  displacement  potentials  xa  are  defined  by  Eq.  (A. 11) 
in  Appendix  A,  and  the  normalization  factor  N . is 


(4tt /e  ) (Jl+m)  ! 

,2  _ m 

Jim  (24+1)  ( J? — m)  ! 


(5.4) 


with  e = 1 and  e = 2,  m ± 0. 
o m 


The  multipole  coefficients  A^  (w)  characterizing  the 
earthquake  source  are  then  calculated  by  taking  the  Fourier 
transform 


<“> 


<k„*> 


jfU 


(Rft) 


, . -iwt 
(R,t)  e dt 


(5.5) 
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The  numerical  calculation  of  the  multipole  coefficients 
in  the  time  domain  expressed  by  Eq.  (5.3)  is  performed  first 
by  the  series  of  general  purpose  computer  program  MULTEES 
described  in  Appendix  G of  Bache,  el:  al.  [1975a]  , The  3-D 
Lagrangian  finite  difference  code  (TRES)  was  carried  into 
the  small  displacement  elastic  region  and  linked  to  MULTEES. 
The  saved  quantities  were  time  histories  for  the  four  poten- 
tials x (a  = 1,  2,  3,  4)  defined  at  the  centers  of  the  eight 
rectangular  cells  surrounding  a total  of  319  monitoring  sta- 
tions. This  set  of  monitoring  stations  was  distributed  over 
the  first  octant  at  grid  points  in  the  neighborhood  of  a 
spherical  surface  of  radius  R = 1.5  km  with  0°  <_  6 <_  90°  and 
0°  < <0  < 90°.  A 3-D  spatial  interpolation  routine  is  applied 
to  derive  values  of  the  potentials  at  a pre-selected  set  of 
points  on  the  spherical  surface.  Once  this  interpolation  is 
accomplished,  the  double  numerical  integration  over  the  two 
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angular  variables  9 and  0 is  performed  as  indicated  in  Eq. 
(5.3).  The  Fourier  transform  given  in  Eq.  (5.5)  is  then 
computed,  yielding  the  multipole  coefficients  which  com- 
pletely specify  the  equivalent  elastic  source. 


5.4  THE  EQUIVALENT  ELASTIC  SOURCE  FOR  THE  BILATERAL  RUPTURE 

The  symmetry  properties  of  the  stick-slip  earthquake 
for  the  configuration  illustrated  in  Fig.  5.1  are  developed 
in  Appendix  A.  These  relations  were  applied  to  the  double 
surface  integrations  defined  by  Eq.  (5.3)  to  derive  the  set 
of  nonvanishing  multipole  coefficients  for  the  bilateral 
rupture  model  (see  Appendix  B) . A summary  of  these  coeffi- 
cients is  given  in  Table  B.l.  The  notation  introduced  above 
for  the  coefficients  is  related  to  that  used  in  earlier 
reports  by  A<“>  5 and  B<£>  e a‘“>.  Only  even  order 

terms  are  found  to  exist  - 0,  2,  4,  6,  . ..);  all  terms 
of  odd  order  are  identically  zero.  There  is  a single  nonzero 
monopole  term  Av^.  There  are  five  quadrupole  terms,  which 
specify  the  double-couple  portion  of  the  equivalent  elastic 
source,  nine  octupole  terms,  etc.  The  spectral  ranges  where 
such  terms  must  be  included  to  ensure  convergence  of  the 
multipolar  series  will  be  a subject  of  later  discussions. 

In  order  to  determine  the  general  nature  of  the 
equivalent  elastic  source  representation  of  the  bilateral 
rupture,  ths  spectral  displacement  potentials  in  the  source 
coordinate  system  can  be  expressed  explicitly  in  terms  of 
the  nonzero  multipole  coefficients  as 
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(R,0,4>,w)  = A(1)(w)  h(2)  (k  R) 

0 0 0 s 

+ [a^  (w)  ?°  (cos  0) 
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(R,6,<fr,w)  = B ^ (to)  P1  (cos  0)  sin  <p  h(2^  (k  R) 
T 2 1 2 2 P 


B^  (<*))  Pl  (cos  6)  sin  <p 

[41  4 

(oj)  P3  (cos  9)  sin  3$  x (k  R) 

4 3 4 4 P 


(5.9) 


where  k = w/v  and  k = w/v  . 

s s p p 

The  radiation  field  is  therefore  made  up  of  monopole 
(Z  = 0)  , quadrupole  (*-  = 2),  octupole  (Z  = 4)  , and  higher 
order  even  multipole  terms.  The  monopole  term  corresponds 
to  a rotation  about  the  x-axis  that  is  the  same  at  every 

point  on  the  surface  of  any  sphere  in  the  elastic  region. 

(1) 

The  time  history  of  the  multipole  coefficient  A (R,t)  indi- 

o o 

cates  that  the  propagating  rupture  produces  primarily  a 
positive  rotation  about  the  x-axis  (in  a positive  rotation  the 
y-axis  rotates  towards  the  z-axis  as  in  the  "right-hand  rule") 
that  is  preceded  and  followed  by  a negative  rotation  about 
the  same  axis.  There  is  no  permanent  offset  so  that  the  net 
effect  of  the  monopole  term  is  an  oscillating  rotation  which 
returns  the  medium  to  its  original  state  at  long  times.  The 
presence  of  the  monopole  term  is  somewhat  unusual,  since 
even  rather  complex  analytical  earthquake  source  models  such 
as  Archambeau’s  [1964]  do  not  generate  such  a term.  Mono- 
pole radiation  is  spherically  symmetric  as,  for  example,  the 
monopole  or  reduced  displacement  potential  representation  of 
an  explosion.  While  perhaps  it  is  not  intuitively  obvious, 
the  excitation  of  S-wave  monopole  radiation  is  a reasonable 
consequence  of  the  model  formulation.  In  the  model  considered 
the  prestress  is  assumed  to  be  pure  shear,  with  only  a single 
non vanishing  component  SZy.  Since  the  fault  plane  is  the  XY 
plane,  a rotational  displacement  potential  with  an  X component 
is  consistent  with  the  prestress  condition. 
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A dominantly  quadrupole  equivalent  elastic  source 
is  expected  for  the  faul-  model  considered.  The  perturba- 
tion by  the  even  multipoles  of  higher  order  is  associated 
with  the  finite  rupture  propagation  velocity.  The  absence 
of  odd  multipoles  is  found  to  be  associated  with  the  bi- 
lateral nature  of  the  rupture. 

The  radiation  pattern  is  specified  by  the  Legendre 
functions  P™  (cos  6)  and  the  trigonometric  functions  cos  m<J> 
and  sin  m<J>.  The  radial  decay  enters  via  the  functions 
h (2)  ^ f which  are  spherical  Hankel  functions  of  the  second 
kind,  while  the  frequency  dependence  of  the  source  is  speci- 
fied  by  the  multipole  coefficients  (u)  . 

The  multipolar  expansion  for  the  equivalent  elastic 
source  is  non-separable.  As  a consequence,  discussions  of 
the  frequency  behavior,  spatial  behavior,  and  the  relative 
contribution  of  the  various  multipolar  terms  cannot  be 
treated  independently.  Conclusions  about  the  nature  of  the 
source  must  then  be  derived  from  displacement  spectra  and 
radiation  patterns  and,  therefore,  must  consider  the  special 
conditions  (spectral  range,  spatial  location  and  series 
truncation)  pertinent  to  the  result  under  consideration. 

The  asymptotic  behavior  of  the  equivalent  elastic 
source  is  described  in  Appendix  C.  In  the  low-frequency 
limit  when  w « 1,  the  far-field  radiation  is  expected  to  be 
dominated  by  the  quadrupole  field,  which  depends  on  the  first 
power  of  a),  namely, 

|xa(R,9,p,w) I * w (5< 

for  k R >>  1 and  u>  <<  1.  The  monopole  and  higher-order  terns 
in  the  potential  spectrum  behave  as  w2,  oj3,  w5,  etc.  The 
asymptotic  behavior  of  the  far-field  displacement  spectrum 
in  the  low-frequency  limit  will  be  independent  of  frequency, 
approaching  a DC  spectral  level. 
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In  the  high-frequency  limit  the  far-field  displacement 
spectrum  is  expected  to  behave  as 


ui  (R,  0 , P/to)  | ^ w 


(5.11 


for  to  >>  1 with  q <_  3.  The  high  frequency  behavior  is  more 
complex,  since  multipole  fields  of  all  orders  contribute  to 
the  radiation.  One  objective  in  analyzing  spectral  amplitudes 
is  to  define  the  frequency  range  where  the  multipole  series 
converges  rapidly,  so  that  retaining  only  the  first  few  terms 
is  sufficient  to  describe  the  radiation  field  at  all  frequencies 
of  interest  at  teleseismic  distances.  High-frequency  disper- 
sion and  errors  due  to  the  discrete  grid  used  in  the  3-D 
finite  difference  calculation  introduce  a high-frequency  cutoff 
(Appendix  D) . The  existence  of  such  a cutoff  limits  the  high- 
frequency  information  that  can  be  extracted  from  the  numerical 
source  calculation.  The  cutoff  frequencies  for  P and  S waves 
for  the  3-D  stress  wave  calculation  are  estimated  to  be 


fc(P)  * 5 Hz 


(5.12 


f (S)  » 3 Hz, 
c 


These  high  frequency  limits  will  not  affect  the  radiation  of 
interest  teleseismically,  where  frequencies  of  interest  are 
generally  below  1-3  Hz. 


FAR-FIELD  DISPLACEMENT  SPECTRA 


In  this  section  the  character  of  the  far-field  dis- 
placement spectra  for  the  finite  difference  bilateral  earth- 
quake calculation  is  discussed.  The  far-field  displacement 
field  is  composed  of  the  radial  or  P-wave  displacement  ur 
and  the  two  components  of  the  S-wave,  Uq  and  u^.  Expres- 
sions for  these  displacements  as  functions  of  R,  8,  <j>,  w are 
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given  by  Eqs.  (C.9  - C.14)  of  Appendix  C.  From  these  ex- 
pressions we  see  that  the  R dependence  is  separable  but 
that  the  dependence  on  the  radiation  pattern  coordinates 
(6,  <j>)  is  not.  Spectra  are  therefore  dependent  on  the  por- 
tion of  the  radiation  pattern  sampled. 

The  coordinate  system  used  for  specifying  position 
with  respect  to  the  fault  is  shown  in  Fig.  5.1.  There  Core, 
for  this  discussion  9 is  the  colatitude  and  is  measured 
from  the  z-axis  while  4)  is  the  azimuth  measured  from  the 
x-axis.  One  usually  thinks  of  the  z-axis  as  being  normal  to 
the  free  surface,  hence  the  nomenclature  for  0 and  4>. 
However,  this  places  the  fault  at  an  unusual  and  not  very 
interesting  orientation  so  this  interpretation  of  the  co- 
ordinate system  is  not  pursued.  Rather,  the  0,  4>  locations 
are  thought  of  with  respect  to  the  source  coordinate  system 
of  Fig.  5.1.  Minster  [1973]  has  given  a computationally  con- 
venient scheme  for  obtaining  the  multipole  coefficients  in 
a rotated  coordinate  system  so  the  fault  studied  here  can  be 
represented  at  any  desired  orientation  with  respect  to  a 
free  surface. 

To  portray  the  far-field  displacement  spectra  radiated 
by  the  FDEQ  it  is  necessary  to  show  radiation  patterns  for 
specific  frequency  components  and  spectra  at  specific  loca- 
tions on  the  radiation  patterns.  These  are  shown  in  Figs. 

5.2  - 5.4.  Figures  5.2  and  5.3  show  a fairly  comprehensive 
sampling  of  radiation  patterns  while  the  three  spectral  plots 
are  typical.  In  the  radiation  pattern  plots  only  the  u0  com- 
ponent of  the  S-wave  is  shown  for  simplicity.  Taken  together, 
these  plots  give  a detailed  understanding  of  the  far-field 
displacement  spectra  radiated  by  the  FDEQ.  In  the  radiation 
patterns  illustrated  in  Figs.  5.2  and  5.3  the  inner  graph 
represents  the  P-wave  pattern.  Compressional  (+)  lobes  and 
rarefaction  (-)  lobes  are  indicated  by  the  plus  and  minus 
signs.  In  the  S-wave  patterns  the  direction  of  transverse 
particle  motion  is  also  indicated  by  arrows. 
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In  Fig.  5.2  a number  of  radiation  patterns  are  shown 
on  the  plane  specified  by  <f>  = 90°;  that  is,  the  ZY  plane. 

In  this  figure  9 = 0°  is  the  +Z-axis,  9 *=  90r  is  the  +Y 
axis,  etc.  The  fault  break  is  along  the  Y-axis.  Recalling 
that  the  displacement  spectra  are  inversely  proportional  to 
R,  the  displacements  in  the  figures  are  scaled  to  1.5  kilo- 
meters, the  elastic  radius.  The  radiation  patterns  are  for 
three  frequencies  spanning  the  tele seismic  band  and  show 
the  influence  of  the  octupole  term.  In  the  top  three  plots 
only  the  monopole  and  quadrupole  terms  were  retained  in 
computing  the  displacements  =0,  2).  In  the  bottom  three 
plots  the  octupole  term  (SL  = 0,  2,  4)  was  included.  From 
Fig.  5.2  we  deduce  the  following: 

• At  0.5  Hz  the  source  is  essentially  a pure 
double-couple.  The  deviations  from  double- 
couple  symmetry  are  quite  small. 

• At  1.0  Hz  the  contribution  of  the  monopole  and 
quadrupole  terms  becomes  evident  in  the  S-wave 
pattern.  At  this  frequency  the  monopole  causes 
a noticeable  asymmetry  in  the  pattern  which  is 
enhanced  by  the  addition  of  the  octupole.  For 
the  more  complete  solution  the  vertically 
radiated  "+"  Ug  waves  are  nearly  1.5  times  as 
large  as  those  radiated  from  the  ends  of  the 
fault.  For  a pure  double-couple  source  the  two 
would,  of  course,  be  equal. 

• At  3.0  Hz  the  monopole  and  octupole  are  larger 
yet  compared  to  the  quadrupole  source. 

• Retaining  terms  through  the  octupole  (£  = 0,  2,  4) 
is  found  to  be  sufficient  for  frequencies  to 

3.0  Hz  or  a little  longer.  Beyond  this, 
higher-order  terms  begin  to  be  felt. 


Ill 
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• Even  for  this  small  earthquake,  terms  other 
than  the  quadrupole  contain  significant  in- 
formation for  far-field  radiation  at  fre- 
quencies observed  teleseismically . 

The  patterns  of  Fig.  5.2  are  all  for  radiation  in  the 
ZY  plane.  A series  of  patterns  at  different  slices  are  shown 
in  Fig.  5.3.  The  patterns  are  all  with  l = 0,  2,  A and  at 
1.0  Hz.  Three  azimuths  (0  = 30° , 60°,  90°)  and  three  colatitudes 
(0  = 90°,  120° , 150°)  are  shown.  These  are  shown  to  indicate 
that  the  conclusions  drawn  above  are  not  influenced  by  the 
particular  azimuth  chosen. 

In  Figs.  5.4  are  shown  spectra  for  all  three  components 
of  displacement  at  representative  locations  on  the  radiation 
pattern.  The  plotted  spectra  are  for  the  azimuths  <J>  = 30°, 

60°,  90°  and  the  colatitudes  0 = 90°,  120°,  150°.  For  these 
spectra,  the  following  comments  are  appropriate: 

• Only  urf  is  radiated  at  0 = 90°;  that  is, 
in  the  XY  plane. 

• The  S-wave  radiation  is  entirely  uQ  at 
<p  = 90°. 

• The  spectral  shape  is  quite  simple  for  all  fre- 
quencies, being  flat  at  low  frequency  and  rolling 
off  smoothly  at  high  frequencies. 

• Inequality  between  the  0 = 120°  and  0 = 150° 
components  of  u„  and  uQ  is  a measure  of 
the  asymmetry  introduced  by  the  monopola  and 
octupole  terms. 

• At  low  frequency  the  source  is  nearly  pure 
quadrupole  or  double-couple. 

• The  spectral  shape  is  essentially  unchanged 
for  all  <p  and  0 shown.  The  spectra  merely 
scale  according  to  their  ratios  at  a single 
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frequency.  That  is,  these  representative 
spectra  together  with  radiation  patterns 
(Figs.  5.2  - 5.3)  completely  describe  the 
displacement  spectra  from  the  FDEQ. 
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Two  other  features  of  the  spectra  are  of  interest,  the 
P/S  amplitude  ratio  and  the  corner  frequency.  For  the  P/S 
ratio  the  radiation  pattern  effect  is  deleted  by  considering 
only  the  peak  values  on  the  respective  P and  S patterns. 

Then  for  this  study  we  define 


(u  (R,  6 = 45°,  * = 90°,  to)) 

P/S  (to)  = 


(5.13) 


(uQ  (R,  e = 0°,  4>  = 90°,  to)) 


At  low  frequencies,  i.e.,  in  the  flat  portion  of  the 
spectra,  this  ratio  is  0.22.  At  1.0  Hz  the  P/S  ratio  de- 
creases to  0.20  while  P/S  (3.0  Hz)  = 0.15.  Therefore,  the 
P/S  ratio  is  constant  at  low  frequencies  and  decreases  at 
higher  frequencies. 

Corner  frequency  is  a somewhat  arbitrary  quantity. 

The  usual  technique  is  to  draw  asymptotes  to  the  high  and  low 
frequency  portions  of  the  spectrum  and  identify  their  inter- 
section with  the  corner  frequency.  Following  this  procedure 
with  the  spectra  in  Fig.  5.4,  we  find  the  corner  frequencies 
to  be  weakly  dependent  on  azimuth  but  to  be  in  the  range 
1.3  - 1.7  Hz  for  both  P and  S waves.  Blandford  [1975]  sug- 
gests that  corner  frequency  be  defined  as  that  frequency  at 
which  the  displacement  spectrum  is  half  its  low  frequency 
limit.  Following  this  procedure,  we  find  the  corner  fre- 
quencies to  be  nearly  independent  of  azimuth  and  to  have  the 
values 


(P) 


= 1.50  Hz 


(£) 


(5.14) 


= 1.55  Hz 
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In  either  case,  the  FDEQ  does  not  exhibit  corner  frequencies 
that  are  substantially  different  for  P and  S waves.  This 
result  is  at  variance  from  the  observation  by  a number  of 
authors  that  the  corner  frequencies  for  P are  higher  than 
for  S. 

The  corner  frequency  is  often  used  to  estimate  the 
earthquake  source  dimension  D using  the  relationship 


(5.15) 


where  Vg  is  the  S velocity.  For  example,  Brune  [1970]  gives 

C„  c = 0.37.  Wyss  and  Shaney  [1975]  recently  determined 
tr  f o 

CD  c for  two  deep  earthquakes  to  be 

C fO 


Cp  = 0.41  - 0.45  , 
Cs  « 0.30  - 0.36  , 


(5.16) 


which  are  consistent  with  Brune' s results. 

For  the  FDEQ  we  have  Vg  = 3.4  km/sec  while  the  fault 

zone  had  a length  and  width  of  1.0  km  and  0.6  km.  If  we  take 

f ■ 1.5  Hz,  C„  . * U.37,  we  find  from  (5.15)  that  D s 0.84  km, 
o * , “ 

a value  which  is  quite  consistent  with  the  fault  dimensions  of 
the  FDEQ. 
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APPEND TX  A 

SYMMETRY  PROPERTIES  OF  THE  STICK-SLIP  RUPTURE  SOURCE  MODEL 


A.l  SYMMETRY  PROPERTIES  OF  THE  PARTICLE  DISPLACEMENT 
FIELL 

The  symmetry  properties  of  the  particle  displacement 
field  are  a consequence  of  the  specific  set  of  boundary  con- 
ditions chosen  for  the  stresses  at  the  fault  plane.  This 
set  of  boundary  conditions  is  required  to  complete  the  defini- 
tion of  the  stick-slip  rupture  model.  In  order  to  check  the 
validity  of  the  boundary  specifications  incorporated  in  the 
3-D  finite  difference  stress  wave  code  (TREL)  , a full-grid 
computation  was  performed  for  a limited  number  (11)  of  time 
cycles.  The  fault  center  was  located  at  the  origin  of  a 
rectangular  coordinate  system,  identified  by  the  cell  indices 

(26.26.26) .  A detailed  examination  was  made  of  the  particle 
displacements  (u,  v,  w) , and  velocities  (u,  v,  w)  at  the  grid 
points  defining  the  eight  cells  surrounding  the  fault  center. 

A comparison  of  the  displacement  and  velocity  components  at 
the  grid  point  (27,27,27)  in  the  first  octant  with  the  points 

(25.27.27) ,  (25,25,27),  (27,25,27),  (27,27,25),  (25,27,25), 

(l 

(25,25,25)  and  (27,25,25)  representing  points  in  the  2nd,  3rd, 
4th,  5th,  Sth,  7th  and  8th  octants,  respectively,  yields  a 
set  of  symmetry  relations  and  boundary  conditions  for  the 
particle  displacement  and  velocity  fields.  The  fault  plane 
lies  in  the  XY  plane  and  rupture  propagates  along  the  Y-axis 
(see  Fig.  5.1). 

Reflection  across  the  YZ  plane  from  the  1st  to  the  2nd 
octant  yields 

u (-x,y , z , t) 
v(-x,y,z, t) 
w ( -x , y , z , t ) 

I I K 


= -u(x,y,z,t)  , 

= +v(x,y,z,t)  , (A.l) 

= +w(x,y,z, t)  , 
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and  in  the  XZ  Plane  <*-«.  the  bounty  condition  is 


u(0,y,z,t)  = u(0,y,z.t)  = 0 . 


(A. 2) 


the  XZ  Plane  fro.  the  2nd  to  the  3rd  octant 
Reflection  across  the  XZ  plan 

yields 


u(-x,-y,Z,t)  = +u(x,y,z,t>  , 

v(-x,-y,z,t)  = +v(x,y.z,t)  , 

w(-x,-y,z,t>  ■=  -„(x,y,z,t)  . 


(A. 3) 


wl-x,-yr ^ 

a the  M plane  from  the  4th  to  the  1st  octant 
Reflection  across  the  XZ  pian 

yields 


u (x,-y »z»t)  = -u(x,y»Zft)  , 

v(x,-y,z,t)  = +v  (x,y , zrt)  , 
W(x,-y,z,t)  = -w (xfy»Zrt)  , 


(A. 4) 


and  in  the  XZ  plane  (y-0)  the  boundary  condrtrons  are 

u(x,0,z,t>  * i(x,0,z,t)  = 0 , (A.5) 

w(x,0,z,t)  ■=  i(x,0,z,t)  - 0 . 

the  XX  Plane  from  the  1st  to  the  5th  octant 
Reflection  across  the  XX  plane 

yields 


u(x,y,-z,t)  = -u(x,y,z,t)  , 
v(x,y,-z,t>  ■=  -v(x,y,z,t)  . 
w(x,y,-z,t)  = +w(x,y,z,t)  . 


(A. 6) 


ties  are  found  relating  displacements  between 
Similar  properta  ^ ^ and  8th  and  4th  octants.  The 

the  6th  and  2n  , are 

boundary  conditions  at  the  XX  pi 

u(x,y,0,t)  = u(x,y,0,t)  = 0 , (n..7) 

V(x,y,0,t)  = vfx.y.O.t)  = o . 
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On  the  Z-axis  (x=y=0) , 
u(0,0,z,t)  = 0 . 
w(0, 0,z,t)  = 0 . 

Since  the  pressure  is  a zone-centered  variable,  values 
of  P at  the  center  of  each  of  the  eight  cells  surrounding  the 
fault  center  at  grid  indices  (26,26,26)  were  examined.  The 
following  symmetry  properties  for  the  pressure  were  derived: 

Octant 

2 P(-x,y,z,t)  = +P(x,y,z,t)  , 

3 P(-x,-y,z,t)  = -P (x,y , z , t)  , 

4 P(x,-y,z,t)  = -P(x,y,z,t)  , 

(A. 8) 

5 P (x,y,-z,t)  = -P (x,y, z , t)  , 

6 P(-x,y,-z,t)  = -P (x,y,z, t)  , 

7 P(-x,-y,-z,t)  = +P(x,y,z,t)  , 

8 P(x,-y,-z,y)  = +P(x,y,z,t)  . 

The  boundary  conditions  observed  for  the  pressure  are 


P(x,0,z,t)  =0  XZ  plane  , 

P(x,y,0,t)  =0  XY  plane  . 


(A. 9) 


The  symmetry  properties  of  the  rupture  model  source  are 
given  for  the  particle  displacement  field  in  Table  A.l.  The 
symbols  u,  v,  and  w represent  the  X,  Y,  and  Z components  of 
che  displacement  field.  Values  for  the  displacement  components 
in  each  octant  are  related  to  the  corresponding  values  in  the 
first  octant.  Consequently,  particle  displacements,  velocities, 
and  potentials  at  all  points  of  interest  on  the  surface  of  the 
sphere  in  the  elastic  region  can  be  derived  from  the  time 
histories  generated  by  the  numerical  source  calculation  for 
the  first  octant. 
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TABLE  A.l.  SYMMETRY  PROPERTIES  OF  THE  DISPLACEMENT  FIELD 

(RECTANGULAR  COORDINATES)  FOR  THE  RUPTURE  MODEL 

(u,  v,  w = X,  Y,  and  Z component  of  particle  displacement) 


Octant 

2 u(-x,y,z)  = -u (x,y,z) 
v(-xry,z)  = +v(x,y,z) 
w(-x,y,z)  « +w(x,y,z) 

3 u (-x,-y, z)  = +u(x,y,z) 
v(-x,-y,z)  = +v(x,y,z) 
w(-x,-y,z)  = -w(x,y,z) 

4 u(x,-y,z)  = -u(x,y,z) 
v(x,-y,z)  = +v  (x , y , z ) 
w(x,-y,z)  = -w(x,y,z) 

5 u (x,y,-z)  = -u(x,y,z) 
v(x,y,-z)  = -v(x,y,z) 
w(x,y,-z)  = +w(x,y,z) 

6 u (-z,y,-z)  = + u (x,y  ,.z) 
v(-z,y,-z)  = “ v (x , y , z ) 
w(-z,y,-z)  = + w(x,y,z) 

7 u(-x,-y,-z)  = -u(x,y,z) 
v(-x,-y,-z)  = ~v (x,y, z) 
w(-x,-y,-z)  = -w(x,y,z) 

8 u(x,-y,-z)  = +u(x,y,z) 
v'x,-y,-z)  = -v(x,y,z) 
w(x,-y,-z)  = -w(x,y,z) 
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While  the  source  symmetry  properties  are  first  defined 
in  rectangular  coordinates  corresponding  to  the  use  of  a 
rectangular  grid  system  in  the  3-D  finite  difference  source 
code,  in  studies  of  the  far-field  displacement  spectra 
spherical  components  are  of  most  interest.  Symmetry  proper- 
ties for  the  displacement  field  in  spherical  coordinates  de- 
rived from  the  following  coordinate  transformation  are  pre- 
sented in  Table  A. 2. 

ur(r,0,<j>)  = sin0  cos<})  u + sin0  sin<}>  v + cos0  w , 

ue(r,0,<j>)  = cos0  cos<}>  u + cos0  sin<J>  v - sin0  w , (A. 10) 

u^  (r , 0 , <$)  = -sin4>  u + cos^  v . 


A.  2 SYMMETRY  PROPERTIES  OF  THE  DISPLACEMENT  POTENTIALS 


The  (scalar)  dilatation  displacement  potential  and  the 
(vector)  rotation  displacement  potential  are  defined  in  terms 
of  the  displacement  field  u(r,t)  by. 


u , 


X = 


u . 


(A. 11) 


The  symmetry  properties  of  the  potentials  can  therefore  be 
derived  from  these  for  the  displacement  field.  The  results 
of  this  derivation  are  given  in  Table  A. 3 for  x and  the 

4 

rectangular  components  of  x»  Since  the  pressure  is  directly 
proportional  to  the  dilatation 

P = -K  V • u , (A. 12) 

where  K is  the  bulk  modulus,  the  symmetry  properties  derived 
for  xk  via  those  for  u,  v,  and  w are  identical  to  those  ob- 
served for  the  pressure  given  in  Eq.  (A. 8). 
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TABLE  A.  2.  SYMMETRY  PROPERTIES  OF  THE  DISPLACEMENT  FIELD 
(SPHERICAL  COORDINATES)  FOR  THE  RUPTURE  MODEL 

(u  , ua,  u.  = spherical  components  of  particle  displacement) 
r y 9 

Octant 

2 ur  (r,0,ir-$)  - + ur  (r , 6 ,$) 

Ug(r,0,ir-if>)  - + u0  (r , 6 , 9 ) 

% (r,6,ir-9 ) - - u^(r,e,<|>) 

3 ur (r , 6 , ir+9)  - - ur(r,0,$) 

ue (r , e , ir+9 > - - u0 (r,e,<f>) 

(r,e,ir+4>)  - - (r,0,$) 

* ur  (r  ,0  - - ur(r,e,«) 

ug  (r,8 ) ■ - Ug  (r,0,$) 

% (r,6,-$)  - + u^(r,0,$) 

5 ur  (r,ir-6,<i)  - - ur(r,0,«(>) 

ufl  (r,ir-8,*i)  - + u0  (r,0,$ ) 
u^(r,ir-0,0)  - - u^(r,0,$) 

6 ur  (r,ir-0,ir-$)  ■ - ur(r,0,$) 

Ug  (r,ir-0,ir-0)  - + Ug(r.,0,$J 

%(r,w-0fir-0)  - + u^(r,0,$) 

7 ur(r,ir-0,ir+i(i)  - + ur  (r,0^) 

Ug  (r,ir-0 , w+9)  - - Ug  (r,0,$) 
u$  (r,ir-6,ir+0)  - + u^  (r,0 , 9) 

8 ur(r, *-«•-*♦)  + ur(r,0,$) 

u0  (r,»-0,-0)  - - Ug(r,0,$) 

u*  (r,tr-0,-0)  « - u (r,0,$) 
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TABLE  A 3 SYMMETRY  PROPERTIES  OF  THE  POTENTIAL  FUNCTIONS  IN 
TABLE  A.  . rectangular  COORDINATES  FOR  THE  RUPTURE  MODEL 


Octant 


(-x,y,z)  * + XH (x,y,z) 


X (-x,-y,z)  - - x (x,y,z) 

4 4 


x (x,-y,z)  - - X (x,y,z) 

4 4 


X4(x,y,-z)  - - (x,y,z) 


X <-x,y,-z)  - - x (x , y , z ) 

4 4 


x4<-x,-y,-z>  - + xk  (x,y,z) 


X„ (x,-y,-z)  = + X (X,y,z) 

4 4 


Xt  (-x,y,z)  - xt  (x,y,z) 


Xj (~x,-y,z)  - + xt (x,y,z) 


Xi (x,-y,z)  ■ + xi  (x,y,z) 


Xt  (x,y,-z)  * + xt  (x,y,z) 


Xt (-x,y,-z)  ■ + xt  (x,y,z) 


X i (-x,-y — z)  * + xt  (x,y,z) 


Xt  (x,-y,-z)  ■ + xt  Cx,y,z) 


t J. 
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TABLE  A. 3.  (Continued) 

2 

X2 (-x,y,z)  - - Xj (x,y,z) 

3 

X2 (-x,-y-z)  « + x2 (x,y,z) 

4 

X2 (x,-y,z)  - - x2 (x,y, z) 

5 

X2 (x,y,-z)  = + x2 (x, y, z) 

6 

X2 (-x,y,-z)  - - X2 (x,y,z) 

7 

X2 (-x,-y,-z)  * + X2 (x,y,z) 

8 

X2 (x,-y,-z)  - - x2 (x,y, z) 

2 

X^ <-x,y,z)  " “ X2 (x,y,z) 

3 

X$ (-x,-y,z)  - - x2 (x,y,z) 

4 

X $ (x,-y, z)  • + Xs (x,y, z) 

5 

X 2 (x,y,-z)  - - X 2 (x,y, z) 

6 

X 2 (-x,y,-z)  - + x 2 (x,y, z) 

7 

X2 (-x,-y,-z)  - + x2 (x,y, z) 

8 

X (x,-y-,z)  ■ - x (x,y,z) 

R-2792 
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Two  examples  will  illustrate  the  method  employed  above. 
From  the  definition  of  the  dilatation  potential, 

X^(x,y,z,t)  = (x,y,z,t)  + ~ (x,y,z,t) 

+ U (x,y,z, t)  . (A. 13) 

d Z 

Applying  the  properties  in  Table  A.l  for  the  2nd  octant 
yields 


X (-x,y,z,t)  = (-x,y,z,t)  + ^ (-x,y,z,t) 

+ (x,y,z,t)  , 

= + x (x,y,z,t) • 

4 

From  the  definition  of  the  rotation  potential, 


X 


i 


\ [curl  u]x 


X (x,-y,z, t) 

i 


1 / 3w  3v  \ 

2 \ YJP yT  " FP zY)  ' 


= x (x,y,z,t) . 


(A. 14) 


(A. 15) 


The  relationships  for  the  other  octants  and  other  potential 
functions  are  derived  in  a similar  manner. 

As  a consequence  of  the  symmetry  properties  for  the 
potential  functions  the  following  boundary  conditions  were 
observed  in  the  time  histories  generated  by  the  3-D  finite 
difference  source  code: 
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1.  On  the  X-Z  plane  corresponding  to  y=0  or 


X (R,8,  0,t)  = 0 

4 


x (R,e,  o , t)  = o . 

2 

2.  On  the  Y-Z  plane  corresponding  to  x=0  or  4>=90 ‘ 


x (R,0,7T/2,t)  = 0 , 
2 


X3  (Rfe, ir/2, t)  = o , 

On  the  X-Y  (fault)  plane  corresponding  to  z=0  or 


0=90' 


x (R,Tr/2,4>,t)  = 0 , 

4 


x (R,TT/2,4>ft)  = 0 , 
3 


4.  On  the  Z-axis  corresponding  to  x=y=0  or  0*0' 


X (R,0,4>,t)  = 0 , 
2 


X (R  f 0 f ^ f t ) = 0 , 


X (RfOr^ft)  = 0 

4 


A.  3 APPLICATION  OF  SYMMETRY  PROPERTIES  TO  BOUNDARY  STATIONS 

A principal  application  of  the  symmetry  properties  for 
the  potential  functions  is  the  determination  of  the  xa  for 
monitoring  stations  on  boundaries  of  the  first  octant  in 
which  the  3-D  finite  difference  calculation  is  carried  out. 

As  described  in  Bache,  et  al . [1975a]  , in  order  to  perform  the 
double  numerical  integration  defined  in  Eq.  (5.3),  values  of 
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X must  be  available  at  selected  integration  points  on  the 
surface  of  the  sphere  of  fixed  radius  in  the  elastic  region. 
The  3-D  finite  difference  code  generates  time  histories  for 
x at  the  centers  of  the  eight  cells  surrounding  each  of  the 
monitoring  or  "save”  stations  defined  at  selected  grid  points. 
The  notation  for  labeling  these  8 cells  is  according  to  Bache, 
et  al.  [1975a].  A cell  is  identified  by  the  three  indices  of 
the  grid  point  with  the  three  highest  integer  values.  The 
station  (I,J,K)  will  sample  xa's  at  the  centers  of  the  follow' 
ing  eight  cells,  in  the  sequence  indicated: 


1.  1+1,  J+l,  K+l 

2.  I,  J+l,  K+l 

3.  I,  J,  K+l 

4.  1+1,  J,  K+l 

5.  1+1,  J+l,  K 

6.  I,  J+l,  K 

7.  I,  J,  K 

8.  1+1,  J,  K 


(A. 20 


The  MULTEES  program  applies  a 3-D  spatial  interpolation  routine 
in  order  to  calculate  values  of  x^  ah  the  integration  points 
of  interest.  The  set  of  xa  generated  by  the  3-D  source  code 
is  deficient,  whenever  the  monitoring  station  is  located  near 
the  XY,  YZ , or  XZ  planes,  which  form  the  boundaries  of  the 
first  octant  in  which  the  xa  are  located. 

If  a monitoring  station  is  on  the  Y-Z  plane  (x=0)  , <t>-90°) 

y are  available  only  for  cells  1,  4,  5 and  8 surrounding  the 
Aa 

sta  ion.  Apdying  the  symmetry  relations  for  the  2nd  octant 
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*cx(2)  " xaU) 


xam  = xa(4) 


a = 1,4 


X0(6)  - X„(5) 


Xa(7)  = Xa  (8) 


(A. 21) 


X„(2)  = - Xa(l) 


Xa(3)  = - x0(4) 


a - >,3 


Xa(6)  = - X0(5) 


Xa(7)  - - Xa<8> 


where  the  integers  in  paient.heses  refer  to  the  cell  numbers 
defined  according  to  the  sequence  in  Eq.  (A. 20). 

If  the  station  is  on  the  X-z  plane  (y=0,  <{>=0o),  x 
are  available  only  for  cells  1,  2,  5 and  6 surrounding  the 

station.  Applying  the  symmetry  relations  for  the  4th  octant 
yields. 


Xa(3)  = Xa<2> 


Xa<4)  * Xa(l) 


a '=  1,3 


X (7)  = x„(6) 
a a 


Xa(8)  = xa<5) 


(A. 22) 


Xa(3)  = -Xs(2) 


XaU>  = -X0(l) 


a = 2,4 


Xa<7)  = -X0(6) 


Xa(8)  = -X„(5) 
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If  the  station  is  on  the  X-Y  plane  (z-0,  0*90°),  Xa 
are  available  only  for  cells  1,  2,  3 and  4 surrounding  the 
station.  Applying  the  symmetry  relations  for  the  5th  octant 

yields , 

Xa<5>  = XaU> 


X„(6)  = Xa<2> 


Xa(7)  = XaO) 


Xa(8)  » Xa(4) 


Xa(5)  = -Xa(D 


Xa(6)  = -Xa(2) 


X (7)  = -X  (3) 
Aa  a 


a = 1,2 


(A. 23) 


a = 3, 


U 


G 


Q 


Xa(8)  = -Xa(4) 


If  the  station  is  on  the  Z-axis  (x=y=0,  6=0°) , x are 


a 

available  only  for  cells  1 and  5.  Values  of  xa  for  cells  2 
and  6 are  obtained  from  those  for  cells  1 and  5 by  applying 
the  rules  given  in  Eq.  (A. 21)  for  the  Y-Z  plane.  Xa  values 
for  cells  3,  4,  7 and  8 are  next  obtained  from  those  for 
cells  1,  2,  5 and  6 by  applying  the  rules  given  in  Eq.  (A. 22) 
for  the  X-Z  plane.  If  a station  is  on  the  Y-axis  (x=z=0, 
6=90°),  4>=90°),  xa  are  available  for  only  cells  1 and  4. 
Values  for  cells  2 and  3 are  obtained  by  first  applying  the 
rules  of  the  Y-Z  plane  and  values  for  cells  5,  6,  7 and  8 
by  next  applying  rules  for  the  X-Y  plane.  If  a station  is  on 
the  X-axis  (y=z=0,  6=90°,  $=0°),  xa  are  available  only  for 
cells  1 and  2.  Values  for  cells  3 and  4 are  obtained  by  first 
applying  the  rules  for  the  X-Z  plane  and  values  for  cells  5, 
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6,  7 and  0 by  next  applying  the  rules  for  the  X-Y  plane. 

In  short,  no  special  rules  are  required  for  stations  on 
either  the  X-,  Y-,  or  Z-axis  if  the  rules  for  the  Y-Z,  X-Z 
and  X-Y  planes,  when  applicable,  are  applied  in  sequence. 

The  FORTRAN  subroutine  INTER3  called  by  the  MULTEES 
program  M3  first  applies  the  rules  def.ned  above  for  any 
station  located  at  a boundary  of  the  first  octant.  The  3-D 
spatial  interpolation  routine  is  then  applied  using  the 
"raw"  xa  at  the  eight  surrounding  cells  and  the  corres- 
ponding shape  functions  defined  in  Section  G.6  of  Bache, 
et  al.  [1975a].  By  means  of  the  above  defined  symmetry 
properties  the  MULTEES  programs  are  able  to  perform  the 
double  numerical  integration  indicated  in  the  definition  of 
the  multipole  coefficients  in  the  time  domain  Eq.  (5.3), 
when  linked  to  a 3-D  finite  difference  calculation  that 
generates  source  variables  only  throughout  a region  of  space 
confined  to  the  first  octant.  Treatment  of  the  3-D  rupture 
source  can  therefore  be  restricted  to  a single  octant  or 
one-eighth  of  the  full  grid  that  would  otherwise  be  required 


u 


o 


The  equivalent  elastic  source  representation  for  any 
arbitrary  seismic  source  is  defined  by  a series  of  multi- 
pole  coefficients  (Eq.  (5.1)).  Calculation  of  these  multi- 
pole coefficients  requires  a numerical  double  integration  over 
the  surface  of  a sphere  of  radius  R in  the  elastic  region. 
The  symmetry  properties  of  the  stick-slip  earthquake  for 
the  configuration  illustrated  in  Fig.  5.1  have  been  described 
above  in  Appendix  A.  In  this  appendix,  these  symmetry  pro- 
perties will  be  applied  to  derive  the  nonvanishing  multipole 
coefficients  for  this  bilateral  rupture  source  model.  As  a 
practical  result  from  the  application  of  these  symmetry 
properties,  the  required  double  integrations  will  be  reduced 
to  the  spherical  surface  subtending  only  a single  octant. 
While  the  nonzero  multipole  coefficients  derive'*  in  this 
manner  define  the  equivalent  elastic  source  representation  of 
the  rupture  model  for  a specific  configuration  in  the  source 
coordinate  system,  the  coefficients  for  any  other  configura- 
tion can  be  derived  from  these  by  performing  the  appropriate 
series  of  transformations  (rotations)  to  the  source  coordi- 
nate system. 

The  expression  for  the  multipole  coefficients  in  the 
time  domain  (Eq.  (5.3))  provides  a convenient  starting  point 
for  the  application  of  symmetry  properties.  The  integration 
in  0 can  first  be  reduced  to 
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N!m|  I ' 


tt/2 


0 *^0 


x (R,0,<t>,t)  + (-) l+m  xft 

a a 


x P™  (cos  0)  sin  9 <30] 

X 


cos  m<f>j 
sin  m<j>| 


d 4) , 


(B.l) 


since 

cos  (tt-0)  = - cos  0 , sin  (tt—  0 ) = sin  0,  and 

p™  (cos  (tt-O))  = '(-)  ^+m  P™  (cos  0)  • 

If  the  potential  xu  is  symmetric  or  antisymmetric  with  respect 
to  reflection  through  the  fault  (X-Y)  plane,  that  is,  if 

xa(R/TT-9,4.,t)  - + xa(R,e,$,t) 


where 

0 <_  0 <_  ti/2  and  0 £ <f>  ^ tt/2  , 
then  the  0 integration  in  Eq.  (B.l)  reduces  to 
ti/2 

f [i  + (-)£+m]  xa(R^^^)  plJ  (cos  e>  sin0  de  • (D,2) 

^ 0 

For  symmetric  potential  functions,  the  following  multipole 
coefficients  vanish 

A?ms  = 0 when  (£+m)  is  odd;  (B.3) 


t 
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similarly,  for  antisymmetric  potential  functions, 


(a) 


A'~'  (R, t)  = 0 when  (2.+m)  is  even. 

Jims 


(B.4) 


The  integration  over  the  aximuthal  angle  <i>  can  also 
be  reduced  to  integrals  over  the  first  octant, 


2tt 


/ 


( cos  m<J> ) 

X^(R,e,<j>,t)  d<J)  , 

( sin  mcf>  J 


a 


tt/2 


= j"  |xa(4>)  + cos  nnr  Xa  ('n-'t1)  + cos  miT  XaU^) 


\ ( cos  m4> 

± X„  (“4>)  J d(f>  ' 

a ) ( sin  ir4>  j 


(B.  5) 


where 


V I 
o 

< ir/2 

/ Xa  ) = xa  (R/°,4>,t) 

cos 

in 

= cos  miT  cos  m<{>  , 

sin 

m 

(tt+<|>) 

= + cos  mir  sin  me})  , 

cos 

m 

(-<f>) 

= cos  m<f>  , 

sin 

m 

(-$> 

= - sin  m4>  • 

(B.6) 


In  the  last  and  third  from  last  relations,  the  positive  signs 
correspond  to  integrals  with  the  factor  cos  m<t>  and  the  negative 
signs  to  those  with  sin  m<{>.  When  m is  even,  the  <J>  integral 
above  can  be  expressed  as, 
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[X0«*>  ± Xft(^“4»)  + Xa(^+4>)  + Xa (-<*») 


cos  m<t» 


sin  m<t» ) 


<a<t>  ; 


(B.  7 ) 


similarly,  when  m is  odd. 


J txa(^)  + Xadr-*)  “ Xa(^)±Xa(^) 


cos  m<J>( 
sin  m<t» ! 


(B.  8 ) 


B.l  SCALAR  DILATATION  POTENTIAL  (x  ) 

4 


For  the  scalar  potential  x the  symmetry  properties 

4 


are  (Appendix  A) : 


x (R,Tr-e,4>,t)  = - x (R,e,$,t)  , 

4 4 


x (r, e ,Tr-ct> , t)  = + x (R,e,<i>,t)  , 

4 4 


x (R,e,Tr+(j>,t)  = - x (R,e,4>,t)  , 

4 4 


X(  (R,e,  -<*>,t)  = - X (R/0,4>,t)  . 

4 4 


(B.  9) 


Since  x is  antisymmetric  in  0, 

4 


= 0 when  (£+m)  is  even, 


(B. 10) 
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Application  of  the  symmetry  properties  expressed  by 
Eq.  (B. 9)  in  the  $ integral,  Eq.  (B.7),  shows  that 


A;4  (R, t)  = 0,  when  m is  even, 

£ms 


A(4)  (R, t ) = 0,  when  m is  odd 

£.ms 

and  s = 0. 


(B. 11) 


Consequently,  the  only  nonvanishing  multipole  coefficients 
for  a = 4 are  expressed  as 


n/2  tt/2 


411  (R-t) 


x (R,e,d>,t) 

4 


x p™  (cos  0)  sin6  d0  sin  m<J>d(}> , 

X/ 


(B. 12) 


where  (&+m)  is  odd,  m is  odd  and  s - 1,  in  short, 


s = 1 


Z — 2,  4,  6,  8, 


m = 1,  3,  ..  (£-1). 

The  nonvanishing  multipole  coefficients  for  a = 4 are 

.(4)  .(4)  (4)  .(4)  (4)  (4) 

A211'  A411'  A431 ' A611 ' A631'  A651'  GtC‘ 


B.  2 VECTOR  ROTATION  POTENTIAL  (x) 


For  the  first  component  of  the  vector  potential,  x r 
the  symmetry  properties  are: 


1 


R-2792 


x (R,  Tr-e , 4>  t t)  = + X (R,6,4>,t), 
1 1 


x (R,e,Tr-<t>,t)  = + x (R,0/<f/t)  , 

i » 


x (R,e#Tr+(j»ft)  = + x (R,e,<f,t)  , 


x (R,0,-<J>,t)  = + x (R,6/<f>,t)  . 

1 1 


Application  of  the  symmetry  properties  in  the  6 and  <j> 
integrals  shows  that 


(B, 13) 


tR-t> 


= 0,  when  U+m)  is  odd, 


Alms  <R't) 


= 0,  when  m is  even  and  s = 1, 


(B. 14) 


A,!^  (R,t)  = 0,  when  :n  is  odd. 

£ms 

Consequently,  the  only  nonvanishing  multipole  coefficients 
for  u = 1 can  be  expressed  as 


tt/2  tt/2 


Atal(R't)=^-  f f X,  (R,e,*,t)  P™  (cos 

Nitm  */0  J0 


6) 


x sin0  d0  cos  m<{>  d4> 


(B. 15) 


where  (£+m)  is  even,  m is  even,  and  s = 0;  in  short, 


s = 0 , 


£ = 0 , 2,  4,  6,  ..., 


m = 0 , 2 , . . . SL  . 


The  only  nonvanishing  terms  for  a = 1 are; 


A(1)  A(1)  A (1)  \U)  AU)  A(1)  etc 

000'  200'  220'  400'  420'  440' 
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For  the  second  component  of  the  vector  potential,  \2> 
the  symmetry  properties  are: 


x (R,  Tt—0  , <j> , t)  = + x (R,e,<{),t)  , 
2 2 


x (R, e,TT-4>,t)  = - x , 

2 2 


x ( R , 0 , tt+<|>  , t ) = + X (R,0»4>ft)  , 
2 2 


x (R,e#-4»ft)  = - x (R#e,<j»,t)  . 

2 2 


From  these  relations,  it  follows  that 


(2) 

l£ms 


(B. 16) 


(R,t) 

= o, 

when 

( &+ m) 

is  odd, 

(R,  t ) 

= 0, 

when 

m is 

even,  s = 0, 

(B. 17 ) 

(R,t) 

= 0, 

when 

m is 

odd. 

(2) 

>S 

(2) 

l£ms 


Consequently,  the  only  nonvanishing  multipole  coefficients 
for  a = 2 can  be  expressed  as 


n/2  tt/2 


411  (R-t)  - 


--  f \ 

n?  Jn  A 


x (R,e,<J>,t)  P™  (cos  0) 


£m  * 0 
x sin0  d0  sin  d0, 

where  (£+m)  is  even,  m is  even,  s = 1;  in  short, 
s = 1 , 

H = 2,  4,  6,  8,  ...  , 
m = 2,  4,  ..  • 

The  only  nonvanishing  terms  for  a = 2 are: 
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A 


(2) 

221' 


(2 

421 


(2 

441 


etc. 


For  the  third  component  of  the  vector  potential,  x 
the  required  symmetry  properties  are 

X (R,TT-e,<|>,  t)  = - x (R,0,<J>,t)  • 

3 3 

x (R,  9 , tt— 4> , t)  = - X (R,0,4>»t)  , 

3 ‘ 3 

X3  (R,0  ,TT  + (J>ft)  = - X 3 (R,  e , , t ) , 
x (R,0,-<{>,t)  = + x (R,0,4>»t)  . 

3 3 

Application  of  these  symmetry  properties  in  the  0 and  <J> 
integrals  shows  that 


(3) 

A£ms  ^R,t^  = when  (^+m)  is  even, 

(3) 

A.  (R,t)  = 0,  when  m is  even, 

£ms 

(3) 

A n (R,t)  = 0,  when  m is  odd  and  s = 1. 
urns 

Consequently,  the  only  nonvanishing  multipole  coefficients 
for  a = 3 can  be  expressed 


tt/2  tt/2 

(R't!  - rr  / / VR'0'*'t>  (»s  e> 

‘‘Jim  0 0 

x sin0  d0  cos  m<J)  d4> , 

where  (£+m)  is  odd,  m is  odd,  and  s = 0;  in  short, 
s = 0 , 

£ = 2 , 4,  6,  ...  , 

m = 1,  3 , . . (£-1)  . 
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u 


u 


o 


Q 


The  only  nonvanishing  terms  for  a = 3 are 


(3 
210' 


a(3)  a(3)  a(3)  (3) 

"430'  A610'  A630 ' A650' 


etc. 


A summary  of  the  nonvanishing  multipole  coefficients 
for  the  rupture  model  vr  th  a prestress  S^y  ^ 0 is  presented 
in  Table  B.l  for  conven-jnt  reference.  The  familiar  nota- 
tion used  in  earlier  reports  (e.g.,  Bache,  ct  al ^ 11974]) 
for  these  coefficients  is  related  to  that  introduced  above 

as. 


A 


(a) 

&mo 


B 


(a) 

Hm 


A 


(a) 

£mi 


(B. 22) 


The  second  column  in  this  table  labeled  " Urns  gives  a unique 

index  which  has  a one— to— one  correspondence  with  the  triple 

set  of  indices  (£,m,s).  Summations  over  l,  m,  and  s are 

equivalent  to  a single  summation  over  the  "£ms"  index,  which 

identifies  all  possible  multipole  terms;  for  example,  all 

terms  such  as  Ag01^  = > which  are  identically  zero,  are 

x,  o 1 * o 

excluded  in  a summation  over  "S-ms". 

The  existence  of  a nonvanishing  S-wave  monopole  term, 
A{1),  should  be  noted.  This  term  is  related  to  a rotation 
about  the  X-axis  that  is  in  the  same  direction  for  all  points 
on  the  surface  of  any  sphere  in  the  elastic  region.  Since 
the  rupture  is  not  allowed  to  induce  a net  rigid  body  rotation 
of  the  medium,  this  term  must  vanish  at  late  times. 

There  are  five  nonvanishing  quadrupole  terms 


and  B 


(4) 
21  * 


These  terms  specify  the  double-couple  characterizing  the  source 
and  are  expected  to  be  the  dominant  terms  for  long  periods. 
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APPENDIX  C 

ASYMPTOTIC  BEHAVIOR  OF  THE  EQUIVALENT  ELASTIC  SOURCE 


The  main  interest  in  analytically  evaluating  asymptotic 
forms  for  the  equivalent  elastic  source  is  twofold.  First, 
these  forms  provide  a check  on  the  numerical  calculations  and 
the  resulting  spectral  plots  and  radiation  patterns.  In 
addition,  the  character  of  the  relative  excitation  of  the 
various  multipoles  representing  the  source  will  illustrate 
how  the  3-D  finite  difference  earthquake  simulation  differs 
from  other  earthquake  models.  The  asymptotic  behavior  of 
the  potential  spectra  will  be  described  first,  followed  by 
a discussion  of  displacement  spectra. 

C.l  ASYMPTOTIC  BEHAVIOR  OF  THE  POTENTIAL  SPECTRA 

C.1.1  Low-Frequency  Behavior 

The  low-frequency  limit  corresponds  to  the  condition 

a)  <<  1.  In  the  series  representation  of  the  potential  spectra 

there  are  two  types  of  frequency  dependent  factors,  A'  (to) 

and  h^  ' (kar) , whose  low-frequency  behavior  will  first  be 

examined  separately.  Numerical  results  from  calculations  of 

the  multipole  coefficients  for  the  finite  difference  earth- 

- 

quake  source  show  that  the  A^^  display  the  following  be- 
havior at  low  frequencies  (even  up  to  1.0  Hz): 


I A(1)  ( u) ) I ■v  u) 3 , w <<  1 

0 0 


I ^Ims (w)  I * ^ ' A > 2,  u>  « 1 , 


(C.l) 


where  again 


A(a) 

xrao 
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The  radial  spreading  of  the  potential  spectra  is  con- 
trolled by  the  Hankel  functions.  Since  the  frequency  enters 
these  functions  only  in  the  argument 


k R = 
a 


a)  R/v 
' a 


f 


one  must  take  into  account  the  hypocentral  distance  R of 
the  observer.  For  the  far-field  radiation  the  appropriate 
asymptotic  form  of  the  Hankel  function  is  given  as 

. £+1 

hl2)  (kcR)  ^ L ET ' <V  » (C'2) 

or 


hj [2)  (kaR)  ^ al 


{kaR  >>  1). 


(C.  3) 


In  the  low-frequency  limit  the  "far-field”  monopole 
portion  of  the  potential  will  behave  as 

|A(1)  (oi)  h(2)  (k  R)  | 'o  oo2  ; 
oo  o a 

all  higher  order  terms  (fl,  = 2,  4,  6,  ...)  will  behave  as 

iA^(u))  h|2)  (kaR)  | 'V'  ^ , i > 2 . (C.4) 

In  the  low-frequency  limit  when  w <<  1,  the  far-field  radiation 
is  expected  to  be  dominated  by  the  quadrupole  field,  which  goes 
as  the  first  power  of  w.  In  summary,  the  asymptotic  behavior 
of  the  "far-field"  potential  spectra  in  the  low-frequency 
limit  is 


ISL  (R/8,4>,w)  l^oi, 

vJt 


for 


k R >>  1 and  w <<  1 . 
a 


(C.  5) 
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While  our  interest  here  is  confined  primarily  to  the 
far  field,  it  can  be  shown  that  the  low-frequency  quadrupole 
field  is  also  expected  to  be  dominant  in  the  "near-field"  zone. 
This  limit  corresponds  to  the  double  condition 

a)  <<  1 and  kaR  <<  1 , 


where  the  appropriate  asymptotic  form  of  the  Hankel  function 
as  the  argument  approaches  0 is  given  by  the  leading  term. 


h<2) <k0R> 


U21)  1 

2lii  (kaR)  4+1 


f 


or 

h^2)(kaR)  a,  of*-1  , (kaR  <<  1)  . (C.6) 

Combining  the  results  given  in  Eq.  (C.l)  and  Eq.  (C.6) 
shows  that,  in  the  low  frequency  "near-field"  limit,  the  monopole 
field  will  behave  asymptotically  as  w2,  while  the  quadrupole 
and  higher  (even)  order  terms  behave  as  u>  . As  a recult, 
the  potential  spectra  in  the  low  frequency  "near-field"  limit 
is  expected  to  behave  as 

|xa(R,0,4>/W)  | J , 


for 

u <<  1 and  k^R  <<  1 . (C.7) 

C.l. 2 High-Frequency  Behavior 

The  high-frequency  limit,  which  corresponds  to  the  con- 
dition on  >>  1,  is  not  as  easy  to  evaluate  due  to  limitations 
on  the  high-frequency  information  generated  by  the  numerical 
source  calculation.  The  basis  for  the  existence  of  a high- 
frequency  cutoff  i;3  discussed  in  Appendix  D.  This  cutoff  in 

147 


\ 


■ - — 


— 


d 


R-2792 


the  finite  difference  calculation,  which  is  estimated  to  be 
about  3.0  - 5.0  Hz,  prevents  our  deriving  any  definitive 
conclusions  about  the  behavior  of  the  multipole  coefficients 
in  the  high-frequency  limit. 

( ct ) 

The  numerical  results  from  calculating  A£mg(w)  up 
to  5 Hz,  however,  does  show  that  the  quadrupole  terms  decay 
at  a rate  given  approximately  as  to  Since  the  Hankel 
function  in  the  far-field  contributes  a term  of  to  \ the 

field  due  to  the  quadrupole  part  of  the  source  might  be  ex- 

-2 

pected  to  behave  as  to  in  the  high-frequency  limit  (to  >>  1)  . 

On  the  other  hand,  the  octupole  terms  were  found  to  behave 
differently  in  the  range  of  1.0  to  5.0  Hz  — not  decaying  as 
rapidly  as  the  quadrupole  terms;  consequently,  no  estimate 
could  be  made  of  the  high-frequency  limit  for  terms  of 
higher  order  than  l = 2. 

In  summary,  multipole  fields  of  all  orders  contribute 
to  the  high-frequency  radiation.  Convergence  of  the  multi- 
pole series  may  not  be  very  rapid  in  the  middle  and  high 
frequency  range;  a large  number  of  higher-order  terms  would 
have  to  be  calculated  in  order  to  derive  the  full  radiation 
field  for  any  given  frequency  beyond  say  5 Hz.  However, 
limitations  at  high  frequencies  will  not  introduce  any 
practical  difficulties  for  those  lower  frequencies,  generally 
below  1 - 3 Hz,  that  are  of  chief  interest  in  teleseismic 
studies. 

C . 2 ASYMPTOTIC  BEHAVIOR  OF  THE  FAR-FIELD  DISPLACEMENT  SPECTRA 

The  displacement  spectrum  is  given  in  terms  of  the  poten- 
tial spectrum  by 


u (r,  to)  = - 


— V x (r,  to) 
kl  “ - 


+ — V X x (r,  to) 

k2  ~ 

KS 


(C.  8) 
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where 


is  the  dilatational  potential  and 


is  the 


rotational  potential.  The  far-field  spherical  (polar)  com- 
ponents of  the  displacement  spectrum  are  according  to 
Minster  [1974] 


[ur  (r,  w)]p 


— h x.  <r'  ' 


k i 3r  “ 


r*e  w)1f  ' “7  lsin4>  3F'  ' cos4>  arH  ' 

ks 


2 

= "7  cos0  cos4>  97^ 


+ cos6  sin<$>  -^r2-  - sin4> 


(C.  9) 


In  the  far-field  the  radial  displacement  results  only  from 
P-wave  terms  involving  X4»  and  the  transverse  displacement 
(uQ  and  u.)  is  only  excited  by  S-wave  terms  involving  x» 

t)  <p  — 

Such  a convenient  separation  of  the  P and  S wave  into  purely 
radial  and  purely  transverse  terms  is  not  possible  in  the 
near-field,  since  each  displacement  component  will  involve 
both  x and  x as  given  in  Eq.  (C.8).  Only  displacement 
radiation  patterns  and  spectral  amplitudes  in  the  far-field 
will  be  discussed  below. 

In  the  far-field  approximation  of  the  radial  deriva- 
tives of  the  potentials  appearing  in  Eq.  (C.9),  only  a single 
term  survives,  irrespective  of  multipole  order,  namely 


h<2)  <v» 


. l 

1 e 


-ikaR 


(C.10) 


T* 


El  !o 


Consequently,  in  the  far-field  the  radial  dependence  of  both 
the  longitudinal  and  transverse  displacement  spectra  can  be 
separated  out  of  the  expressions  given  above,  leaving  a 
series  in  which  each  term  is  the  product  only  of  a multipole 
coefficient  and  an  associated  Legendre  function..  For 
example,  the  radial  component  of  the  P-wave  in  the  far-field 
is 


i / 


i I 


U 


u 


u 


o 


o 


. 


o 


1 — 


u^P)  (R,  0 , 0 ,uj)  = 


1_  e 

u R 


-ik  R 
P 


i-1 


i2,  (u>)  P™(cos6)  sin  m<J> 


(C.ll) 


1=2 , 4 , . . m=l 
(odd) 


The  0 component  of  the  S-wave  in  the  far-field  is 


-ikgR 


S<S)  (R,  9 ,4> » w)  = 


KR 


/ °°  x, 

\ A^  (w)  P™(cos6)  cos  m<J)  sin<J> 

v £=0,2,4  m=0 

(even) 


00  x 

y (w)  P™(cos0)  sin  m<j>  cos<}>  f . 


1=2 ,4,6  m=2 

(even) 


) 


(C. 12) 


Expanding  the  far-field  displacement  spectrum  yields, 


■V  i=rr-:=:  ^aa 


I 
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L>  ! 


U*P)  (R, 6 , 4> , uj)  = + g-g---- 

v 2 K 
*P 


* I ^21^  ^ P2  (cos0)  sin<|>  - (to)  P^(cos0)  sin 4) 


S(4) 


- (to)  P*(cos0)  sin34>  + ...  , 


(C. 13 ) 


G<s)  (R,6#4)  ,to) 


-ik  R 
s 


x AqQ  (to)  sin <))  - A^g  (w)  P®(cos0)  sin4> 


\ n 2 


- Aj2  (id)  P^cosS)  cos24>  sin4> 


+ (oi)  P*  (cos0)  sin24>  cos4> 


+ A(1)  (to)  P°  {cos0 ) sin4> 
40  q 


+ A^  (w)  P*(cos0)  cos24>  sin4> 


•**  ( 1 ] . 

+ A^  (to)  Pj  <cos0)  cos44>  sin4> 


(w)  P*(cos0)  sin24>  cos4) 


4 ^ (<J°)  P4(cos6)  sin44>  cost})  + 


(C. 14) 


Analyzing  either  the  P-  or  S-wave  displacement  spectrum 
in  the  far  field  allows  one  to  derive  its  asymptotic  behavior. 


The  low  frequency  asymptotic  behavior  of  the  multipole  co- 

£ 

efficients  has  been  seen  to  be  to  for  Jt  > 2.  The  leading 


[■ 
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term  at  low  frequencies  is  the  quadrupole,  so  that  the  multi- 
pole coefficient  behaves  as  id1.  The  leading  factor  involv- 

- 2 

ing  the  wave  number  k behaves  as  o>  , The  net  result  is 
that  in  the  far-field  the  displacement  spectrum  should  be 
"flat"  in  the  low  frequency  limit.  Only  at  high  frequencies 
when  the  effects  of  multipoles  beyond  1 = 2 become  significant 
will  this  flat  behavior  be  perturbed. 

The  high  frequency  asymptotic  behavior  of  the  displace- 
ment spectrum  is  more  difficult  to  analyze.  While  informa- 
tion on  the  high  frequency  behavior  of  the  multipole  coeffi- 
cients is  limited  in  the  finite  difference  calculation  (see 
Appendix  D) , most  quadrupole  terms  were  estimated  to  behave 
approximately  as  <d-1.  Higher  order  terms  are  not  expected 
to  damp  as  rapidly  as  the  quadrupole  terras,  so  that  the  net 

result  might  be  that  the  asymptotic  behavior  of  the  displace- 

— 2 

ment  spectrum  is  <d  for  <d  >>  1.  Consequently,  the  displace- 
ment spectrum  is  expected  to  behave  as 

|u.  (r,oj)  | 'v*  (D-q  for  u >>  1 , (C.15) 

i ~ 


with  q _<  3. 

The  main  interest  in  analytically  evaluating  asymp- 
totic forms  for  the  spectra  is  that  they  will  provide  a check 
on  the  numerical  calculations  to  be  presented  below.  The 
slope  of  the  amplitude  spectra  is  expected  to  be  azimuthally 
dependent  because  of  interferences  between  the  multipole 
fields  of  different  orders. 
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HIGH-FREQUENCY  APPROXIMATIONS  DUE  TO  THE  DISCRETE 
GRID  IN  THE  THREE-DIMENSIONAL  FINITE 
DIFFERENCE  CALCULATION 

There  are  several  effects  associated  with  the  use  of 
a discrete  grid  which  must  be  taken  into  account  when  analyzing 
the  results  of  numerical  finite  difference  calculations.  A 
practical  result  of  these  effects  is  the  introduction  of  a 
high-frequency  cutoff.  Useful  information  from  the  3-D  earth- 
quake simulation  will  be  limited  to  frequencies  below  this 
cutoff.  Before  proceeding  to  derive  an  estimate  of  the  high- 
frequency  cutoff,  which  is  the  objective  of  this  discussion, 
two  effects  will  be  sketched  to  illustrate  how  approximations 
originate  in  the  finite  difference  method.  High-frequency 
attenuation  due  to  the  introduction  of  artificial  damping 
w*  11  be  discussed  first.  The  second  illustration  will  show 
how  dispersion  effects  arise  from  the  discretization  of  the 
applicable  continuum  mechanics  equations  in  the  3-D  finite 
difference  calculation.  Finally,  the  last  subsection  will 
provide  some  estimates  for  a high-frequency  cutoff. 

D.l  HIGH-FREQUENCY  ATTENUATION  DUE  TO  ARTIFICIAL  DAMPING 

IN  THE  THREE-DIMENSIONAL  EARTHQUAKE  CALCULATION 

In  finite  difference  stress  wave  calculations  it  is 
often  desirable  to  include  viscosity  in  the  material  proper- 
ties to  maintain  smooth  profiles  at  wave  fronts.  This  vis- 
cosity is  often  called  artificial  viscosity.  However,  it 
is  useful  to  realize  that  finite  difference  codes  generally 
do  not  propagate  waves  in  purely  elastic  media,  but  only  in 
materials  that  are  to  some  extent  viscoelastic. 

In  the  three-dimensional  finite  difference  earthquake 
calculation,  the  "artificial  viscosity"  terms  cause  the 
material  to  behave  as  a Voight  solid  model  of  viscoelasticity. 
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The  effect  of  the  viscoelasticity  on  a propagating  wave  is 
shown  in  Ewing,  Jardetsky  and  Press  [1957]  , Chapter  5.  A 
propagating  compressional  wave  is  attenuated  by  exp  [-tx] 

where 


2 _ 0) 

2a2 (1  + a2w2 ) 


1 + a2w2  1/2  - 1 » 


(D.l) 


where  a is  the  compressional  wave  velocity  while  a is  a 
selected  constant,  characterizing  the  material  viscoelasticity. 

In  the  finite  difference  code  a is  selected  to  be 
q At  where  At  is  the  time  step  and  q is  a number  which 
is  empirically  chosen  to  maintain  stability.  It  is,  of 
course,  directly  related  to  the  amount  of  damping. 


When  a2w2  <<  1,  (D.l)  can  be  expanded  in  power  series 


to  give 

t=^(a)2  +0  (a  2 u)  ** ) . 

For  a finite  difference  earthquake  calculation 
a = 5.7  x 10 5 cm/sec  and  a * 10" s sec.  In  this  case, 


(D.2) 


X ~ 8.8  x 10 


-10  w2 


(D.3) 


The  elastic  radius  in  the  earthquake  calculation  was  at  1.5 
km.  Since  the  attenuation  due  to  linear  artificial  viscosity 
is  exp  [-  tx] , its  effect  on  the  displacement  spectrum  can 
be  estimated.  Tabulating  versus  frequency  we  find 


f (Hz) 


exp[-T*1.5  km] 


0.995 


10.0 


0.98 

0.95 

0.92 

0.88 

0.59 
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Therefore,  the  attenuation  due  to  artificial  viscosity  has 
minor  effect  (at  this  distance)  for  frequencies  less  than 
5 Hz. 


D.2  DISCRETIZATION  ERRORS 

The  difference  equations  used  in  finite  difference 
codes  are  discrete  analogs  of  differential  equations.  How- 
ever, the  solutions  of  the  difference  equations  are  not 
identical  to  the  solutions  of  the  analogous  differential 
equations.  We  must  be  certain  that  the  errors  introduced 
by  the  discretization  are  not  influencing  the  solutions  to 
a substantial  degree. 

The  governing  equations  for  one-dimensional,  linear 
elastic  wave  propagation  may  be  written: 

pv  = a'  , 

( 

o = Kv ' , 


(D.4) 


where  v is  velocity,  a is  stress,  K = A + 2p , the  dot 
signifies  time  differentiation  and  the  prime  signifies  spatial 
dif ferentation.  The  finite  difference  equations  in  most 
Lagrangian  finite  difference  codes  are  the  (multi-dimensional) 
finite  difference  analogs  of  the  differential  equations  (D.4). 

We  can  write  (D.4)  in  finite  difference  form  and  study 
the  properties  of  a propagating  wave  compared  to  the  same  wave 
in  a continuous  medium.  The  wave  in  the  discrete  system  turns 
out  to  be  dispersed.  We  should  be  certain  that  the  frequency 
components  of  interest  are  in  the  region  where  this  disper- 
sion effect  is  weak. 

The  finite  difference  analog  of  (D.4)  may  be  written 


p m+1  _ m _ 1 m m 

j Vj  Ax  aj+l/2  " ° j-1/2 


(D.  5) 


1 m+1  _m  \ m+1  m+1 

At  j +1/2  j+1/2  lx  j+1  Vj 
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If  solutions  of  the  form  exp[i (kx-wt) ] are  assumed, 
substitution  into  (D.5)  gives  two  homogeneous  simultaneous 
equations  in  the  unknowns  v and  a.  Requiring  the  deter- 
minant to  vanish  gives  the  dispersion  relation  between  to 
and  k. 


This  is 


p -2itoAt  -i'oAt  , , Ke 

rr-r  e - 2e  + 1 + 


-itoAt 


ikAx 


-ikAx  n 
- e = 0. 


(D.6) 


Note  that  At  and  Ax  are  related  by  the  Courant  stability 


condition 


nAx  = a At, 


(D.7) 


where  n < 1. 


If  it  is  assumed  that  toAt,  kAx  <<  1,  the  exponentials 
in  (D.6)  can  be  expanded  in  power  series  and  the  result  re- 


duced to 


= k2a2  1 + yj  w2At2  - + 0 (<o4At4) 


+ 0 (k4 Ax4 ) , 


(D.8) 


where  a2  = K/p  has  been  used. 
Then  using  (B.7) 


to2  = k2uz  1 + k2 a2  At2  ( 


-— ) 
12n2  / 


(D.  9) 


In  terms  of  the  phase  velocity  c = to/k, 


c = a 1 + ^ to2 At 2 


(D. 10) 


1 •'  ■■ 


. 
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where  b = (25  - n_2)/12  and  for  0.2  < n < 1,  0<b<2. 


\ 


■ 


For  the  three-dimensional  finite  difference  earthquake 
calculation,  the  pertinent  parameters  were  At  = 0.005  sec 
and  n = 0.285.  Then 


C. 


O 


o 


(5 


c s a [1  + 5.2  x lCT**  f 2]  , (D.ll) 

and  we  see  that  the  effect  is  very  small  for  low  frequencies. 

For  5 Hz,  c ~ 1.01  a. 


D. 3 HIGH-FREQUENCY  CUTOFF  ESTIMATES 

Effects  due  to  the  use  of  a discrete  grid  in  the  finite 
difference  calculation  will  now  be  discussed  from  the  point  of 
view  of  deriving  an  estimate  of  the  high-frequency  cutoff  that 
is  associated  with  the  3-D  earthquake  source  calculation.  Con- 
sider first  the  effect  often  referred  to  as  the  aliasing  of 
the  continuous  time  series.  Due  to  the  digitization  of  the 
stress  wave  signal,  information  is  not  available  at  frequencies 
corresponding  to  wave  numbers  greater  than  the  Nyquist  wave 
number;  this  gives  a cutoff  frequency 


fN  - 


v k.. 
a N 

2n 


a 

2 Ax 


where  kN  - 


(D.12) 


is  the  Nyquist  wavenumber,  Ax  is  a grid  size  (cell  length), 
and  vq  the  body  wave  propagation  velocity. 

Another  problem  encountered  by  short  duration  pulses 
is  that  the  velocity  used  to  advance  tiie  motion  at  time  t to 
t + At  is  dependent  on  the  wavenumber  k for  stress  waves 
propagating  on  a discrete  lattice.  The  origin  of  this  dis- 
persion effect  was  sketched  above.  The  dispersion  becomes 
important  for  10  or  fewer  grid  points  per  wavelength  (Boore 
[1972]).  This  limitation  yields  a cutoff  frequency 


157 


V V 
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The  same  result  is  found,  if  one  considers  that  at  least  10 
grid  points  are  required  to  define  a high  frequency  pulse. 

In  a similar  approach,  if  one  uses  the  rule  that  in  a single 
time  cycle  the  wave  can  travel  at  most  only  one-half  a cell 
dimension,  then 


~ = vAt,  where  At  = cycle  time 
2 


and  if  say  10  grid  points  are  required  to  define  the  shortest 
wavelength  then 


**  ■ ir-=  (it)/<10Ax>  = 2k  • 

m i n 


(D. 14) 


In  the  finite  difference  source  calculations  the  following 
parameters  were  used: 


Ax  = 0.1  km 


At  = 0.005  sec 


vp  = 5.7  km/sec 


vc  = 3.4  km/ sec  . 

u 


Substituting  in  the  above  equations  gives  the  following  fre- 
quency limits  (in  Hz) : 


f N = 28*5 


fN  = 17*5  ' 


V = 5.7 
max 


f =3.4,  and 
max 


f = 10  . 
T 


As  a conservative  limit,  for  this  specific  finite  difference 
source  calculation  the  cutoff  frequencies  for  P-wave  and 


* 
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S— wave  spectra  can  be  taken  as  5.7  and  3.4  Hz,  respectively. 
Consequently,  information  should  be  derived  only  for  fre- 
quencies below,  say 


t - 5 Hz 


and 


3 Hz 


CD. 15) 


